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ABSTRACT 


The  objective  of  the  research  project  was  to  numerically  predict  the  characteristics  of 
evaporating  sprays  in  a  turbulent  shear  layer.  The  calculations  were  performed  such  that  the 
essential  physical  processes  in  spray  evaporation  in  gas  turbine  combustors  are  adequately  realized. 
The  calculations  modeled  parallel  experiments  (ARO  Contract  number  DAAG29-84-K-0165) 
carried  out  in  a  vertically  down-flowing  wind  tunnel  where  a  splitter  plate  separates  two  air  streams 
having  different  velocities.  The  two  air  streams  meet  at  the  end  of  the  splitter  plate  and  give  rise  to  a 
confined  shear  layer.  A  flat  prefilming  airblast  atomizer  is  located  at  the  tip  of  the  splitter  plate.  A 
liquid  spray  formed  at  the  atomizer  tip  evaporates  in  the  turbulent  shear  layer  and  is  convected 
downstream. 

An  Eulerian  description  of  tiie  turbulent  gas  phase  flow  and  a  Lagrangian  formulation  of 
the  droplet  motion  and  heat  and  mass  transfer  were  considered.  Single  phase  two  dimensional 
turbulent  flow  calculations  in  a  duct  were  considered.  A  parabolic  formulation  of  the  governing 

equations  along  with  a  k-e  turbulence  model  was  used  to  obtain  the  computed  flow  fields.  The 
values  of  turbulent  mixing  time  within  the  shear  layer  were  computed  and  compared  with  those 
predicted  by  the  characteristic  time  model.  For  droplets  injected  into  the  turbulent  air  stream,  a 
deterministic  separated  model  was  utilized  for  calculating  their  trajectories,  temperatures  and 
evaporation  rates.  Typically,  a  Rosin-Rammler  size  distribution  function  was  employed  for  the 
droplets  at  the  inlet.  Solutions  were  obtained  for  both  water  and  Jet-A  fuel  as  the  liquid  phase.  The 
calculated  SMD  (Sauter  mean  diameter)  variations  along  the  tunnel  axis  compared  favorably  with 
the  measurements.  To  account  for  the  turbulent  dispersion  of  the  droplets  due  to  the  gas  phase 
turbulence,  a  stochastic  separated  flow  model  was  developed.  In  particular,  the  spread  rate  of  the 
spray  due  to  the  dispersion  effects  were  evaluated. 


3 


TABLE  OF  CONTENTS 

Page 

Abstract .  3 

Table  of  Contents .  4 

Acknowledgements . . .  6 

Nomenclature .  7 

I.  INTRODUCTION  AND  BACKGROUND .  9 

II.  ANALYSIS  OF  TURBULENT  SHEAR  FLOW  IN  A  DUCT 

II.  1.  General  approach .  13 

11.2.  Turbulence  model .  14 

11.3.  Initial  and  boundary  conditions  for  the  mean  flow  quantities .  15 

11.4.  Initial  and  boundary  conditions  for  the  turbulent  flow  quantities..  17 

II.  5.  Numerical  Methods .  20 

II.  6  Results  and  discussion .  21 

IB.  ANALYSIS  OF  DROPLET  BEHAVIOR -DETERMINISTIC  MODEL 

III.  1  General  approach . : .  43 

113.2  Momentum  equations .  44 

111. 3  Heat  and  mass  transfer  equations  for  droplets .  45 

111. 4  Droplet  class  and  distribution  functions .  53 

111. 5  Numerical  Methods .  54 

111. 6  Results  and  discussion .  55 

IV.  ANALYSIS  OF  DROPLET  BEHAVIOR  -  STOCHASTIC  MODEL 

IV.  1  Background  and  general  approach .  60 

IV.2  Mathematical  model .  63 

IV. 3  Modifications  to  the  PTRAK  code .  64 

IV. 4  Results  and  discussion .  68 

V.  SUMMARY  AND  CONCLUSIONS .  73 


4 


VI.  REFERENCES 


75 


APPENDICES 


1.  List  of  all  publications  under  sponsorship  of  ARO 

2.  List  of  all  participating  scientific  personnel 


5 


The  author  wishes  to  thank  Professor  A.  M.  Mellor  for  providing  the  initial  thoughts  and 
encouragements  to  carry  out  the  research  performed.  Thanks  also  go  to  Mr.  K.  V.  Tallio  for  his 
interest  in  the  work  .  A  special  note  of  thanks  to  Messrs.  W.  Lau  and  X.  X.  Cai  for  devoting  many 
hours  on  this  project 

The  support  of  the  Army  Research  Office,  and  in  particular  the  assistance  provided  by  Dr. 
David  M.  Mann  as  Technical  Monitor  is  appreciated. 


NOMENCLATURE 


Drag  coefficient  for  the  droplets 
Specific  heat 

Empirical  constant  for  the  turbulence  model 

Droplet  diameter 
Duct  width 
rate  of  strain 

Droplet  heat  transfer  coefficient 
Turbulent  kinetic  energy 
Mass  transfer  coefficient 
Thermal  conductivity 

Turbulent  length  scale,  characteristic  length 


t 

M 

P 

Pr 

O 

q 

R 

Re 

SMD 

Sc 

T 

t 

U 

V 
u’ 
v' 

o 

W 

X 

y 

Y 


mixing  length 
Molecular  weight 
Pressure 
Prandtl  number 
Droplet  heat  transfer 
Universal  gas  constant 
Reynolds  number 
Sauter  mean  diameter 
Schmidt  number 
Temperature 
Time 

Mean  axial  velocity 
Mean  normal  velocity 
Fluctuating  (axial )  velocity 
Fluctuating  (normal)  velocity 
Droplet  evaporation  rate 
Axial  coordinate 
Normal  coordinate 

Fuel  vapor  concentration  (mole  fraction) 


7 


8 

e 

X 

P 

t 

a 

P 

Q 

rm 

Subscripts 

a 

d 

f 

I. 

s 

w 


Boundary  layer  thickness 

Dissipation  rate  of  turbulence 

Shear  layer  strength,  heat  of  vaporization  of  fuel 

Density 

Shear  stress,  relaxation  time 
Standard  deviation 
Viscosity 
Fuel  flow  rate 

Mass  diffusivity  of  fuel  vapor 


Air 

Droplet 

Fuel 

Liquid  phase 
Surface 

Wall  conditions 


8 


I. 


INTRODUCTION  AND  BACKGROUND 


An  experimental  research  program  was  initiated  at  Drexel  University  under  ARO 
sponsorship  (Grant  Number  DAAG29-84-K-0165)  for  continuing  validation  of  the  characteristic 
time  model  [Mellor,  1980].  The  characteristic  time  model  is  a  simplified  view  of  combustion  where 
each  of  the  physical  processes  which  occurs  during  combustion  is  considered  and  a  physical  time 
which  corresponds  to  that  is  computed.  The  semi-empirical  characteristic  time  models  were 
developed  under  TARADCOM/ARO/MERADCOM  sponsorship  for  prediction  of  the  effects  of  fuel 
properties,  combustor  inlet  conditions  and  liner  geometry  on  performance  of  gas  turbine  engines. 
Models  are  available  for  the  correlation  and  prediction  of  combustion  efficiency,  gaseous  emissions, 
lean  flame  stabilization  and  ignitability  (cold  start  and  altitude  relight)  in  terms  of  characteristic  times 
thought  to  describe  the  fluid  mechanics,  chemistry  and  fuel  droplet  evaporation.  By  analyzing  data 
from  helicopter  and  aircraft  engines  it  was  possible  to  generate  performance  curves  using  the 
characteristic  time  model  which  appear  to  be  'universal',  that  is  independent  of  engine  and  fuel  type. 
This  is  believed  to  occur  because  all  conventional  turbine  combustors  utilize  turbulent  spray 
diffusion  flames  for  their  heat  release.  At  ignition,  the  lean  blowoff  limit  and  other  low  power 
operating  conditions,  fuel  evaporation  also  becomes  important.  For  evaporating  sprays  in  a 

turbulent  shear  layer,  the  only  significant  characteristic  times  are  xsj  (the  residence  time  in  the  shear 
layer)  and  xe^  (the  droplet  evaporation  time). 

In  the  above  program  (Grant  Number  DAAG29-84-K-0165),  specially  designed 
laboratory'  tests  were  performed  in  which  the  turbulent  mixing  processes  were  examined  in  detail  for 
comparison  with  these  models.  The  objectives  were  primarily  to  explore  the  characteristic  time 
model  developed  for  gas  turbine  combustors.The  fuel  evaporation  and  turbulent  mixing  times  where 
the  emphasis  of  the  program  where  each  could  be  varied  independently  in  fundamental  experiments. 
A  vertical  down-flowing  wind  tunnel  (see  Figure  1)  with  a  square  test  section  and  provisions  for 
laser  optics  (for  droplet  size  measurements  by  the  forward  diffraction  technique  and  evaporation 
measurements  by  the  absorption  technique)  and  hot  wire/film  probe  insertion  were  designed  and 

constructed  for  local  measurements  of  xsj  and  Xe^.  Presence  of  a  splitter  plate  which  terminates  in  a 

flat  prefilming  airblast  atomizer  produced  a  shear  layer  and  fuel  spray  within  the  tunnel  .The 
situation  in  the  outer  primary  zone  region  of  a  conventional  gas  turbine  combustor  thus  could  be 
closely  approximated  in  the  test  section. 

The  experimental  configuration  selected,  unlike  combustors,  is  two  dimensional  (as 
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shown  by  r  .casurements),  provides  optical  access,  easy  probe  (hot  wire/film)  access,  a  liquid  spray 
injected  uito  a  shear  layer  and  a  region  of  primary  zone  combustor  flow  that  is  important  to  many 
aspects  of  its  performance.  Further,  the  measurements  are  related  to  fundamental  interactions 
between  turbulent  mixing  and  liquid  evaporation  (without  combustion)  in  a  clean  flow  which  can  be 
analyzed  in  an  orderly  fashion  without  all  of  the  complexities  of  practical  hardware.  The 
cross-sectional  dimensions  of  the  tunnel  are  7.65  cm  by  7.65  cm,  chosen  to  provide  test  section 
mean  velocities  ranging  from  50  to  130  m/s.  The  measurements  were  conducted  near  atmospheric 
pressure,  however  liquid  and  air  preheats  to  400  K  and  450  K  could  be  provided.  Air  to  the 
laboratory  rig  flowed  from  a  Reliance  multi-stage  axial  blower  and  heating  was  provided  by  a 
Chromalox  in-line  electrical  resistance  air-heater.  This  system  was  designed  to  deliver  an  air  mass 
flow  rate  of  1.17  kg/s  at  maximum  blower  output. 


Figure  1.  Overall  system  schematic 


As  noted  earlier,  a  flat  prefilming  airblast  atomizer  (with  a  porous  plug  at  the  end)  was 
located  inside  the  test  section  of  the  wind  tunnel.  Upstream  of  and  adjacent  to  the  injector,  a  splitter 
plate  physically  separated  the  two  settling  chambers,  one  for  flow  in  each  side  of  the  splitter 
plate/injector.  This  arrangement  allowed  different  air  velocities  on  the  two  sides  of  the  injector  and 
establishes  a  shear  layer  emanating  from  the  nozzle  tip  as  shown  in  Figure  2.  The  average  flow 
along  the  side  of  the  porous  plug  is  labeled  as  the  'fuel  side’  (Uafs)  and  the  average  flow  along  the 
other  side  is  labeled  as  'air  side'  (Uaas).  The  difference  between  the  velocities  at  the  two  sides  at 
the  tip  of  the  injector  determines  the  strength  of  the  shear  layer  being  studied.  In  addition  to 
boundary  layers  that  develop  along  the  duct  walls,  sharp  velocity  gradients  are  also  present 
downstream  of  the  splitter  plate  due  to  a  shear  layer  or  a  wake. 


Figure  2.  The  shear  layer  formed  within  the  test  section 


In  the  present  research  project,  detailed  calculations  for  the  evaporating  sprays  were 
performed  to  complement  the  above  experimental  program.  The  experimental  and  numerical  results 


provide  a  basic  simplified  test  bed  for  the  the  validation  of  the  chaaeteristic  time  model.  Finite 
difference  calculations  for  the  turbulent  flow  fields  in  a  confined  shear  layer  were  carried  out. 
Calculations  for  the  trajectories  and  evaporating  characteristics  for  injected  droplets  were  then 
carried  out  using  a  deterministic  separated  flow  model.  In  addition,  the  present  numerical 
calculations  were  extended  to  give  insights  to  the  turbulence-spray  interactions  viz.  the  dispersion 
of  liquid  droplets  due  to  gas  phase  turbulence. 

Two  dimensional  steady  state  predictions  of  the  turbulent  shear  flows  within  a  vertical 
duct  (similar  to  the  experimental  configuration)  were  obtained.  The  predictions  were  used  to  validate 
the  experimental  measurements  in  addition  to  extending  the  range  of  the  parameters  considered  in 
the  experimental  studies.  Analysis  in  the  present  investigations  was  limited  to  dilute  sprays;  typical 
of  gas  turbine  combustors.  The  finite  difference  calculations  of  the  flow  field  and  spray  behavior 
allows  evaluation  of  the  local  characteris  tic  times  of  turbulent  mixing  and  droplet  evaporation  for 
comparison  with  both  semi-empirical  [Mellor,1980]  and  the  experimental  values. 


II.  ANALYSIS  OF  TURBULENT  SHEAR  FLOW  IN  A  DIJCT 
II.l.  General  approach 

Assuming  dilute  sprays,  calculation  of  the  two  dimensional  turbulent  flow  field 
(neglecting  the  effects  of  fuel  injection  and  the  fuel  insertion  device)  was  carried  out  first.  The  ADD 
code  [Anderson  et  al.,  1982J  was  utilized  to  perform  the  calculations  .  The  ADD  code  along  with 
the  PTRAK  and  VAPDIF  codes  are  among  a  series  of  three  computer  codes  developed  for 
predicting  the  distribution  of  liquid  fuel  droplets  and  fuel  vapor  in  premixing  -prevaporizing  fuel -air 
mixing  passages  of  the  direct  injection  type.  For  a  given  flow  prediction  obtained  by  the  ADD  code, 
the  PTRAK  code  predicts  the  injected  fuel  droplet  trajectories  and  evaporation  char*jferistics  using  a 
deterministic  separated  flow  model.  The  VAPDIF  code  (not  used  in  the  present  research  program) 
predicts  the  vapor  diffusion  based  on  the  droplet  evaporation  characteristics  calculated  by  the 
PTRAK  code.  In  applying  this  approach,  it  is  implicitly  assumed  that  the  air  flow  behavior  can 
affect  the  fuel  droplet  behavior,  but  the  fuel  droplet  behavior  does  not  affect  the  air  flow  behavior. 
This  is  justified  when  the  mass  fractions  of  fuel  droplets  and  fuel  vapor  are  small. 

The  ADD  code  was  developed  to  solve  the  internal  flow  weak  interaction  problem  using  a 
forward  marching  numerical  procedure  that  does  not  require  an  interaction  between  the  inviscid  core 
flow  described  by  the  Euler  equations  and  the  boundary  layers  described  by  the  parabolic  boundary 
layer  equations.  A  forward  marching  scheme  is  used  to  solve  the  gas  phase  equations  for  mass, 
momentum  and  energy.  An  orthogonal  coordinate  system  is  constructed  for  the  duct  from  the 
potential  flow  solution  such  that  the  stream  function  forms  the  coordinate  normal  to  the  wall  and  the 
velocity  potential  forms  .the  cwdinate  tangent  to  the  wall.  Since  the  potential  flow  streamlines 
approximate  the  real  streamlines  (in  a  two  dimensional  sense),  the  equations  of  motion  may  be 
greatly  simplified  by  assuming  that  the  velocity  normal  to  the  potential  flow  streamlines  is  small 
compared  to  the  streamwise  velocity.  This  procedure  reduces  the  governing  viscous  flow  equations 
to  a  parabolic  system  of  partial  differential  equations  which  can  be  solved  by  a  forward  marching 
numerical  integration  procedure.  This  procedure  is  attractive  when  a  long  tunnel  needs  to  be 
simulated.  For  the  present  calculations,  all  calculations  were  done  for  a  rectangular  (two 
dimensional)  flow  domain  with  straight  duct  walls.  Hence  the  orthogonal  streamline  coorj/nates 
generated  by  the  ADD  code  were  identical  to  the  Cartesian  (x  -  y )  coordinates. 

The  basic  parabolized  compressible  Navier  Stokes  equations  solved  by  the  ADD  code  are 
given  by  Anderson  and  Edwards  [  198 1 J  and  are  not  repeated  in  this  report. 
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II.2.  Turbulence  model 

For  prediction  of  the  turbulent  shear  flow,  time-averaged  conservation  equations  were 
considered  in  the  analysis.  Various  turbulence  models  have  been  proposed  for  the  closure  of  the 
time  -averaged  equations  in  the  past.  One  disadvantage  of  time  averaging  is  that  some  important 
details  of  turbulent  flows  (e.g  coherent  structures)  are  lost  in  the  predictions.  However,  depending 
on  the  sophistication  of  the  turbulence  model  used,  mixing  characteristics  and  measurable  time 
averaged  quantities  (e.g.  the  turbulence  intensity  and  the  length  scales  )  can  be  predicted  with  the 
time-averaged  conservation  equations  with  reasonable  computing  times.  Some  of  these  quantities 
have  been  used  in  the  calculation  of  characteristic  mixing  times  in  turbulent  shear  layers  [  Mellor, 
1980]. 


The  algebraic  turbulence  models  based  on  Prandtl's  mixing  length  theory  are  va  lid  for 
equilibrium  turbulent  boundary  layer  and  wake  flows.  In  the  algebraic  models,  the  eddy  viscosity  at 
any  point  in  the  flow  and  hence  the  turbulent  Reynolds  stresses  are  a  function  only  of  the  local 
mean  flow  and  not  on  the  flow  history.  One  equation  turbulence  models  such  as  described  by 
Gibson  et  al.  [1970]  and  the  two  equation  models  of  turbulence  such  as  those  formulated  by 
Launder  and  coworkers  [1977]  were  developed  to  treat  nonequilibrium  flows  in  which  the  turbulent 

Reynolds  stress  is  a  function  of  the  flow  history.  A  two  equation  model  ( k  -  e  )  incorporated  into 
the  ADD  code  was  used  for  the  present  predictions.  The  model  was  originally  developed  by  Jones 
and  Launder  [1972],  modified  for  streamline  curvature  effects  by  Launder  et  al.  [1977]  and  then 

modified  by  Chen  [1980]  for  low  Reynolds  numbers.  The  turbulent  viscosity  in  the  k-e  model  is 
given  as 

px  =  C^Pk2/e  (II.  1) 

where  is  an  empirical  constant  .  The  model  contains  other  empirical  constants  in  the  transport 
equation  for  the  dissipation  rate  of  turbulence. 

The  governing  equations  for  the  turbulent  kinetic  energy  (  k )  and  its  rate  of  dissipation 

(e  )  in  the  streamline  orthogonal  coordinate  system  (as  required  for  the  ADD  code  formulation)  are 
given  in  Anderson  and  Edwards  [1981].  For  most  of  the  calculations,  the  empirical  constants  used 

in  the  model  are  as  those  given  by  Chen  [1980],  A  sensitivity  analysis  was  done  few  the  constant  C^ 
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for  better  agreement  of  the  predicted  mean  and  turbulent  gas  flow  quantities  with  the  measured 
values  in  the  parallel  experiments  (see  Section  11.6). 

II.3.  Initial  and  boundary  conditions  for  the  mean  flow  quantities 


A  parabolic  formulation  of  the  governing  equations  along  with  a  k-e  turbulence  model 
was  used  to  obtain  the  computed  flow  fields  in  the  duct.  A  forward  marching  scheme  is  used  to 
numerically  calculate  the  flow  field.  The  splitter  plate  (as  shown  in  Figure  2)  was  not  included  in  the 
calculations.  The  mean  and  turbulent  flow  quantities  at  the  end  of  the  splitter  plate  are  prescribed  as 
the  inlet  boundary  conditions. 

For  the  mean  flow  conditions  at  the  inlet,  the  x  component  (streamwise)  of  the  mean 
velocity  distribution  needs  to  be  prescribed  .  The  y  component  (normal)  of  the  mean  velocity  is  not 
solved  in  the  ADD  code  formulation  (orthogonal  streamline  cordinate)  and  hence  no  inlet 
distribution  is  required  for  it.  Any  arbitrary  distribution  of  the  streamwise  velocity  distribution  can 
be  prescribed.  In  addition  to  mean  velocities,  static  pressure  and  temperature  of  the  incoming  flow 
are  also  prescribed.  For  the  duct  walls,  no-slip  conditions  for  the  velocity  is  used.  The  walls  are 
considered  to  be  insulated.  For  the  calculations  presented,  the  inlet  static  pressure  and  temperature 
of  the  two  streams  (on  the  two  sides  of  the  splitter  plate)  are  considered  to  be  the  same.  Mean  axial 
velocity  profiles  measured  close  to  the  end  of  the  splitter  plate  are  used  as  the  inlet  conditions  for 
obtaining  flow  predictions  downstream.  For  equal  average  velocities  of  two  streams  (Uafs  =  Uaas) 
on  the  two  sides  of  the  splitter  plate,  the  inlet  velocity  profiles  are  e.'  .racterized  by  a  wake  region 

(due  to  the  thickness  of  the  splitter  plate).  For  non-equal  average  velocities  (Uafs  *  Uaas)  a  shear 
layer  is  present  at  the  central  region  in  addition  to  the  boundary  layers  near  the  duct  walls. 

Modifications  to  the  ADD  code 

The  ADD  code  was  written  for  analyzing  (straight  or  curved)  duct  flows  with  boundary 
layers  near  the  duct  walls  and  a  core  flow  region  near  the  center.  To  resolve  the  sharp  velocity 
gradients  near  a  wall,  a  turbulent  boundary  layer  profile  is  calculated  in  the  ADD  code  from  the 
prescribed  inlet  mean  velocity  distribution  [  Anderson  et  al.,  1982],  The  boundary  layer  profile 
utilized  in  ADD  [Coles,  1956]  is  an  empirical  formulation.  The  shape  of  the  turbulent  boundary 
layers  also  used  in  grid  distribution  along  the  normal  direction.  In  general,  this  creates  a  dense  grid 
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distribution  near  the  duct  walls  and  a  coarse  distribution  near  the  central  region  of  the  duct.  For  the 
present  flow  configuration  (a  confined  turbulent  shear  layer  produced  by  a  splitter  plate),  dense 
grids  were  required  at  the  central  region  ,  as  well  as  near  the  duct  wall  regions.  In  order  to  have 
sufficient  number  of  grid  points  near  the  central  region  of  the  duct,  a  grid  stretching  parameter  (a 
built-in  feature  in  the  ADD  code)  was  utilized.  This  was  suitable  for  most  of  the  computations 
performed.  However,  for  very  strong  shear  layers  at  the  central  region  of  the  duct,  the  grid  was 
redistributed  at  the  core  region  to  increase  the  number  of  nodes  at  the  central  region  of  the  duct.  This 
was  done  without  disturbing  the  grid  distribution  near  the  wall . 

A  typical  mesh  (uncorrected)  produced  by  the  ADD  code  (for  a  prescribed  inlet  flow 

distribution  with  a  shear  layer  at  the  central  zone)  is  shown  in  Figure  3  where  13  x  55  (x  x  y)  grid 
is  used.  Figure  4  gives  the  redistributed  computational  mesh  for  the  confined  shear  layer 
calculations. . 


X/D 

figure  3.  Mesh  produced  by  the  ADD  code  from  the  inlet  flow  distribution 
measured  at  the  end  of  the  splitter  plate. 
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Figure  4.  Modified  mesh  for  the  same  initial  flow  distribution  data  as  in  Figure  3. 

Most  of  the  results  presented  in  this  report  were  obtained  with  a  redistributed  mesh  as  the 
inlet  flow  distribution  at  the  inlet  contained  a  shear  layer  or  a  wake  at  the  central  region. 


II.4.  Initial  an 


UamifriaBl 


For  the  two  equation  turbulence  model  used  in  the  study,  distributions  for  k  and  e  are 
required  at  the  exit  of  the  splitter  plate.  These  inlet  values  are  not  known  a  priori,  even  if  the  mean 

axial  flow  distribution  is  specified.  For  the  present  calculations,  the  inlet  values  for  k  and  e  were 
obtained  either  from  the  mixing  layer  theory  [Anderson  and  Edwards,  1981]  or  derived  from 
experimental  measurements  of  the  turbulence  intensities  and  length  scales  close  to  the  end  of  the 

splitter  plate.  At  the  duct  walls,  the  values  of  k  and  E  were  both  set  equal  to  zero  along  the  entire 
length. 
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The  Prandtl  mixing  length  tj  defined  by 

HT  =  Pt2|el  (H.2) 

where  px  is  the  turbulent  viscosity  and  e  is  the  rate  of  strain.  In  the  region  near  the  wall,  the 
mixing  length  is  modified  by  the  van  Driest  factor  (A+)  which  is  given  by 

pwuVmv  -  kY+[1  -exp(-Y+/A+)]  (IL3) 

where  u*  is  the  friction  velocity 

V^w  /pw 

Y+  is  the  universal  distance  to  the  wall  (pw  u*y  /  p^)  and  K  is  the  von  Karman  constant. 

In  the  outer  region  of  the  turbulent  boundary  layer,  including  the  free  stream, 

MT  =  P  i  umax  (H.4) 

where  the  mixing  length  is  given  by 

t  =  X§*  (II.5) 

5*  is  the  displacement  thickness  of  the  turbulent  boundary  layer  and  the  proportionality  constant  x 

is  determined  by  Clauser  [1956]  to  be  equal  to  0.016.  The  initial  conditions  for  k  and  e  are  then 
given  as  (within  the  boundary  layer): 


k  =  ypTIX|/(qi0-5.  pE.  p) 

and 

(II.6) 

e  =  Y  pX  1 22  |  /  ( PE2.  p ) 

(II- 7) 

where  pg  =  p  +  px  and  y  is  an  intermittency  factor  given  by  Klebanoff  [1954]  to  reduce  the 

turbulence  to  zero  at  the  edge  of  the  boundary  layer: 

Y  JU  5.5  (|-)6] 

fa .  s) 

The  boundary  layer  thickness  5  is  taken  as  halfwidth  of  the  duct. 


I  is  the  streamwise  stress  distribution  at  the  inlet  (obtained  from  the  prescribed  or  measured  inlet 


flow  profile.  For  the  present  flow  calculations  E  reduces  to  (due  to  straight  duct ) 

L  (D.9) 

where  U(y)  is  the  mean  axial  flow  distribution  at  the  inlet 

For  a  duct  flow  without  a  shear  layer  or  a  wake  (i.e.  no  splitter  plate  upstream)  ,  the 
central  region  generally  has  a  flat  velocity  profile  (free  stream  condition).  Though  such  cases  were 
not  considered  in  the  present  study,  the  turbulent  kinetic  energy  and  dissipation  rate  can  be  specified 
at  the  inlet  for  such  cases  if  the  turbulence  intensity  is  specified.  Thus  if  I  is  the  turbulence  intensity, 
the  turbulent  kinetic  energy  is  given  as 

k--£uL,  <II10) 

The  free  stream  dissipation  is  taken  as  the  value  e  in  the  boundary  layer  where  k  in  the  boundary 
layer  matches  koo. 


Evaluation  of  k  and  e  at  the  inlet  from  the  measurements  of  the  turbulence  intensities  and  length 

scales  at  tbe.end  of  the  splitter  plats: 

Assuming  isotropic  turbulence,  the  k  and  e  values  at  the  inlet  of  the  computational  domain 
can  also  be  evaluated  from  the  measured  distributions  of  the  turbulence  intensity  (rms  fluctuations) 
and  the  longitudinal  length  scale  Le.  These  vlaues  were  measured  (see  the  final  repert  to  ARO 
Grant  Number  DAAG29-84-K-0165)  close  to  the  splitter  plate  end  at  x/D  =  0.03  (see  Figure  2)  as  a 
function  of  the  y  direction. 

The  turbulent  kinetic  energy  is  defind  as 

k  =  0.5  (u'2  +  v'2  +  w'2  )  01.11) 

where  u'  ,  v'  and  w'  are  the  fluctuating  velocity  componenets  at  the  x,  y  and  z  directions.  For 
isotropic  turbulence. 
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(IU2) 


k  =  2/3  u’2 


From  Jones  and  Launder  [1972],  the  length  sacle  may  be  derived  from  the  definitions  of  k  and  e 
as: 

Lg  =  Cjj3/4  k^/2  /  e  (n.13) 

from  where  the  dissipation  rate  can  be  expressed  as 


e  =  C^3/4  k3/2  /  Le 


01.14) 


Hence  from  the  measured  values  of  turbulent  intesities 


and  Le  along  the  y  direction  at  the  end  of  the  splitter  plate,  the  initial  conditions  for  k  and  e  can  be 
specified. 


II.5.  Numerical  Methods  ’ 


The  relationship  between  (lg  and  the  mean  flow  are  specified  through  the  parabolized 

governing  equations  (Anderson  and  Edwards,  1981)  and  are  solved  by  a  forward  marching 
numerical  integration  scheme.  The  equations  are  first  linearized  by  expanding  all  dependent 

variables  in  a  Taylor  series  expansion  in  the  marching  direction  (x)  and  terms  of  O  (Ax^)  are 
dropped.  Finite  difference  equations  are  then  obtained  using  the  two  point  centered  difference 
scheme  of  Keller  [1970, 1968].  The  resulting  matrix  equations  are  block  tridiagonal  and  are  solved 
by  block  factorization  using  the  method  of  Varah  [1972]. 

The  numerical  solution  is  second  order  accurate  in  the  y  direction,  first  order  in  the  x 
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direction  and  linearly  stable.  The  Ax  step  size  is  limited  not  by  linear  stability  conditions  but  by 
required  accuracy  in  the  Taylor  series  expansion  in  x. 


The  numerical  procedure  for  solving  the  governing  equations  for  k  and  e  is  as  follows. 

The  procedure  involves  solving  the  k  and  e  transport  equations  successively  in  an  iterative 
procedure  at  each  stream  wise  station  and  lagging  the  mean  flow  equations  one  step  so  that  the  mean 
flow  variables  are  always  known.  Furthermore,  some  of  the  nonlinear  finite  difference  source  terms 
are  arranged  on  the  central  diagonal  of  the  tridiagonal  matrix  to  strengthen  matrix  stability 
[Anderson  and  Edwards,  1981]. 

II.6  Results  and  discussion 


The  computations  carried  out  for  the  turbulent  shear  flows  within  the  duct  closely  parallel 
the  series  of  experimental  measurements  performed  under  ARO  Grant  No.  DAAG29-84-K-0165. 
Table  1  lists  the  cases  for  which  measurements  were  done  for  the  mean  gas  velocities,  the 
turbulence  intensity  (rms  velocity), 

V? 

and  the  length  scale  Le  at  various  axial  locations. 
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Table  1 

Gas  phase  (experimental )  matrix 


CASE 

Uafs 

(m/sec) 

Uaas 

(m/sec) 

X 

1 

91.4 

91.4 

0.0 

2 

91.4 

18.1 

0.67 

3 

91.4 

0.0 

1.0 

4 

54.8 

27.6 

0.33 

5 

91.4 

46.0 

0.33 

6 

122.0 

61.5 

0.33 

7  ' 

54.8 

54.8 

0.0 

8 

122.0 

122.0 

0.0 

10 

105.0 

105.0 

0.0 

11 

105.0 

85.9 

0.1 

12 

105.0 

70.0 

0.2 

13 

122.0 

122.0 

0.0 

14 

91.4 

60.9 

0.2 

Detailed  measurements  of  mean  and  turbulent  quantities  for  the  above  cases  can  be  found  in  the 
final  report  of  ARO  Grant  No.  DAAG29-84-K-0165.  Uafs  and  Uaas  denote  the  mean  (spatially 
averaged)  gas  velocities  in  the  two  sections  of  the  duct,  separated  by  the  splitter  plate  (see  Figure  2). 
X  is  an  indication  of  the  strength  of  the  shear  layer  defined  as 


X  = 


u 


afs 


The  inlet  conditions  for  the  computations  were  obtained  from  measurements  of  U(y)  at 
x/D  =  0.03.  The  initial  conditions  for  k  and  e  can  be  evaluated  from  the  inlet  mean  velocity 
distribution  using  the  mixing  layer  theory  (as  outlined  in  an  earlier  section)  or  they  can  be  obtained 
from  the  turbulence  intensity  and  length  scale  measurements  performed  at  x/D  =  0.03.  Mean  gas 
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velocities ,  turbulence  intensity  and  length  scale  distributions  obtained  from  the  turbulent  kinetic 
energy  and  dissipation  rate  predictions  were  compared  with  the  experimental  measurements 
performed  at  specific  downstream  locations.  The  computations  were  not  performed  for  all  cases 
listed  in  Table  1  due  to  the  parabolic  formulation  of  the  governing  equations  in  the  ADD  code.  In 
particular,  for  cases  2  and  3,  the  shear  layer  strength  is  high  and  it  is  likely  that  reverse  flow 
situations  may  exist  at  the  tip  of  the  splitter  plate.  Computed  results  for  cases  1,  3,  4,  5,  6, 7  and  8 
are  presented  in  this  report. 

The  calculaions  are  initiated  at  the  exit  plane  of  the  splitter  plate.  A  parametric  study  was 
carried  out  to  investigate  the  effects  of  adjusting  the  measured  length  scale  values  at  x/D  =0.03  on 
the  prediction  of  turbulence  intensities  downstream.  The  measured  length  scale  values  for  case  7 
was  used  for  this  purpose.  The  parametric  study  was  performed  by  using  the  unadjusted  measured 
length  scale  values,  one-half  the  length  scale  values  and  one-quarter  of  the  length  scale  values  (for 
evaluating  the  initial  values  of  dissiapation  rate  of  turbulent  kinetic  energy  e).  Comaprisons  of  the 
predictions  with  experimental  measurements  of  turbulence  intensities  were  made  at  x/D  =  1.  It  was 
observed  that  the  predicted  intensities  were  directly  proportional  to  the  length  scale  sizes  at  inlet. 
The  best  agreement  of  the  predictions  with  the  experiments  were  obtained  when  the  inlet  length 
scale  values  were  1/2  of  the  measured  values  near  the  core  flow  and  1/4  of  the  measured  values  near 
the  shear  layer  (extending  from  approximately  y/D  =  0.15  to  y/D  =  -  0.15).  The  above  adjustment 
can  be  justified  from  the  inherent  limitations  in  the  computational  model  resulting  from  the 
assumption  of  isotropic  turbulence  in  the  duct  flow.  All  results  presented  in  the  section  were 
obtained  by  using  the  measured  inlet  length  scale  distributions,  with  the  above  adjustment. 

A  relationship  between  e  and  the  length  scale  is  provided  by  equation  (11.14): 

e  =  C^3/4  k3/2  / 

The  value  of  is  traditionally  held  constant  and  is  set  equal  to  0.09  in  k  -  e  turbulence  models. 
Sonin  [1983]  reported  that  the  value  of  could  be  calibrated  and  depends  on  the  nature  of  the 
flow  regime  (wake,  free  jet,  boundary  layer  etc.).  A  parametric  study  was  undertaken  to  obtain  a 
value  for  that  gave  best  agreement  between  the  predictions  and  experiments  performed  at 
downstream  locations  of  the  splitter  plate  exit  plane.  The  results  indicated  that  a  value  of  = 
0.571  (corresponding  to  reducing  the  value  of  Le  by  a  factor  of  1/4)  gave  the  best  agreement 
between  experimental  and  computational  profiles  of  turbulence  intensity  in  the  shear  layer  region. 
However,  a  value  of  =  0.045  provided  better  agreement  in  the  wall  boundary  layer.  It  was  also 
observed  that  that  the  turbulence  intensity  predictions  were  insensitive  to  values  of  used  in  the 
core  region  of  the  flow  where  velocity  gradients  are  smaller. 


23 


Figures  5a,  5b  and  5c  show  the  comparison  of  the  predicted  and  measured  mean  axial 
velocities  along  the  duct  width  at  the  inlet,  x/d  =  1  and  2  respectively,  for  case  1.  The  measured 
values  at  x/d  =  0.03  are  used  to  generate  the  inlet  velocity  profiles  for  the  computations.  Similar 
plots  for  turbulence  intensity 


at  the  inlet,  x/D  =1  and  2  are  shown  in  Figures  6a,  6b  and  6c  for  case  1.  Umax  is  the  maximum 
local  velocity  measured  at  the  inlet  (x/d  =  0.03).  Comparison  of  the  measured  and  predicted  mean 
velocities  and  turbulence  intensities  along  the  width  of  the  duct  are  shown  in  Figs.  7a,  7b, 7c  and 
8a,  8b,  8c  for  case  4,  in  Figs.  9a,  9b,  9c  and  10a,  10b,  10c  for  case  5,  in  Figs.  11a,  lib,  11c  and 
12a,  12b,  12c  for  case  6,  in  Figs.  13a,  13b,  13c,  and  14a, 14b, 14c  for  case  7  and  in 
Figs. 15a, 15b, 15c  and  16a,16b,16c  for  case  8.  The  square  symbols  represent  the  computed  results 
while  the  triangles  denote  experimental  measurements 


The  predicted  mean  velocity  profiles  agreed  favorably  with  the  measured  data  in  most 
instances,  Previous  computations  with  uncorrected  grids  (see  Section  II. 3)  resulted  in  poor 
agreement  of  measured  and  predicted  values.  Computational  results  of  the  mean  velocity  profile  at 
the  center  consistently  fell  behind  the  measured  values.  This  is  thought  to  be  due  to  the 
simplifications  made  in  the  parabolic  formulation  of  the  conservation  equations.  The  predicted 
turbulence  intensity 

“ 
u _ 


profiles  were  obtained  from  the  turbulent  kinetic  energy  predictions  using 

.  3  ~ 

k  =  y  u 

The  agreement  of  the  predictions  with  the  measured  values  is  good,  considering  the  inherent 
limitations  of  the  k  -  e  model  employed  in  the  calculations. 
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Figure  5a.  Comparison  of  the  mean  velocity  profiles  at  inlet  for  Case  1 
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Figure  5b.  Comparison  of  the  mean  velocity  profiles  at  x/D  =  1  for  Case  1 
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Figure  6a.  Comparison  of  the  turbulence  intensities  at  inlet  for  Case  1 
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Figure  6b.  Comparison  of  the  turbulence  intensities  at  x/D  =  1  for  Case  1 
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Figure  6c.  Comparison  of  the  turbulence  intensities  at  x/D  =2  for  Case  1 
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Figure  7  a.  Comparison  of  the  mean  velocity  profiles  at  inlet  for  Case  4 


Figure  7b.  Comparison  of  the  mean  velocity  profiles  at  x/D  =  1  for  Case  4 
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Figure  7c.  Comparison  of  the  mean  velocity  profiles  at  x/D  =  2  for  Case  4 
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Figure  8a.  Comparison  of  the  turbulence  intensities  at  inlet  for  Case  4 
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Figure  8b.  Comparison  of  the  turbulence  intensities  at  x/D  =  1  for  Case  4 


Figure  8c.  Comparison  of  the  turbulence  intensities  at  x/D  =  2  for  Case  4 
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Figure  9a.  Comparison  of  the  mean  velocity  profiles  at  inlet  for  Case  5 
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Figure  9b.  Comparison  of  the  mean  velocity  profiles  at  x/D  =  1  for  Case  5 
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Figure  9c.  Comparison  of  the  mean  velocity  profiles  at  x/D  =  2  for  Case  5 
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Figure  10a.  Comparison  of  the  turbulence  intensities  at  inlet  for  Case  5 
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Figure  10b.  Comparison  of  the  turbulence  intensities  at  x/D  =  1  for  Case  5 
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Figure  10c.  Comparison  of  the  turbulence  intensities  at  x/D  =2  for  Case  5 
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Figure  1  la.  Comparison  of  the  mean  velocity  profiles  at  inlet  for  Case  6 
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Figure  1 1  b.  Comparison  of  the  mean  velocity  profiles  at  x/D  =  1  for  Case  6 
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Figure  1  lc.  Comparison  of  the  mean  velocity  profiles  at  x/D  =  2  for  Case  6 
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Figure  12a.  Comparison  of  the  turbulence  intensities  at  inlet  for  Case  6 
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Figure  12b.  Comparison  of  the  turbulence  intensities  at  x/D  =1  for  Case  6 
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Figure  12c.  Comparison  of  the  turbulence  intensities  at  x/D  =  2  for  Case  6 
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Figure  1 3a.  Comparison  of  the  mean  velocity  profiles  at  inlet  for  Case  7 
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Figure  13b.  Comparison  of  the  mean  velocity  profiles  at  x/D  =  1  for  Case  7 
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Figure  13c.  Comparison  of  the  mean  velocity  profiles  at  x/D  =  2  for  Case  7 


Figure  14a.  Comparison  of  the  turbulence  intensities  at  inlet  for  Case  7 
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Figure  I4b.  Comparison  of  the  turbulence  intensities  at  x/D  =1  for  Case  7 


Figure  14c.  Comparison  of  the  turbulence  intensities  at  x/d  =  2  for  Case  7 
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Figure  15a.  Comparison  of  the  mean  velocity  profiles  at  inlet  for  Case  8 
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Figure  15b.  Comparison  of  the  mean  velocity  profiles  at  x/D  =  1  for  Case  8 
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Figure  15c.  Comparison  of  the  mean  velocity  profiles  at  x/D  =  2  for  Case  8 
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Figure  16a.  Comparison  of  the  turbulence  intensities  at  inlet  for  Case  8 


41 


Ur-ns  U>.  trs.:) 


42 


III. 


m.i  General,  approach 

The  behavior  of  the  fuel  spray  in  the  duct  downstream  of  the  splitter  plate  (see  Figure  2) 
is  calculated  using  the  PTRAK  computer  program  [Anderson  et  al.,1982].  The  PTRAK  code  is 
based  on  a  deterministic  separated  flow  model  (mean  gas  flow  variables  are  only  employed  for  the 
droplet  evaporation  characteristics  calculation}  .The  code  was  significantly  modified  so  that  the  fuel 
spray  behavior  can  be  predicted  using  a  stochastic  separated  flow  model  (as  described  in  the  next 
chapter).  The  droplet  dispersion  effects  due  to  gas  phase  turbulence  was  investigated  using  the 
stochastic  separated  flow  model. 

The  present  analyses  assumes  lean  equivalence  ratios,  i.e.,  the  mass  fraction  of  fuel 
(dispersed  phase  )  is  small.  Under  these  conditions,  the  air  flow  behavior  is  solved  first  using  the 
method  described  in  earlier  sections  .  The  pressure  ,  temperature  and  mean  velocity  fields  were  then 
stored  in  data  files  and  used  in  the  fuel  droplet  analysis.  Thus  the  air  flow  behavior  influences  the 
fuel  droplet  behavior,  but  the  fuel  droplet  effects  on  the  gas  flow  are  neglected. 

The  spray  formation  processes  at  the  airblast  atomizer  are  not  addressed  in  the  present 
study.  It  is  assumed  that  the  liquid  layer  on  the  porous  plug  (at  the  end  of  the  splitter  plate) 
disintegrates  within  a  short  distance  and  a  spray  is  formed  consisting  of  droplets  that  can  be 
characterized  by  a  suitable  size  distribution.  For  the  present  calculations,  the  spray  is  assumed  to  be 
formed  at  the  inlet  of  the  computational  domain.  The  spray  is  divided  into  a  number  of  classes.  For 
each  class,  droplet  diameter,  droplet  temperature,  veocity  components  and  position  in  the  flowfield 
are  specified.  These  parameters  may  differ  among  the  various  classes.  The  PTRAK  model 
calculates  the  trajectory,  vaporization  rate  and  temperature  variation  of  each  class  of  droplets.  A 
maximum  of  1250  classes  can  be  considered  in  the  PTRAK  code. 

Droplet  shattering  and  coalescence  were  not  considered  in  the  present  study.  These 
features  are  present  in  the  PTRAK  code  but  are  largely  untested  and  without  validation.  When 
droplets  collide  with  a  duct  wall,  the  predictions  can  be  obtained  using  any  of  the  two  following 
schemes.  Elastic  rebounds  can  be  considered  or  the  heat  transfer  from  the  wall  to  the  droplets  can  be 
considered.  Both  water  (a  pure  compund)  and  Jet-A  fuel  (a  distillate  liquid  )  were  considered  as  the 
dispersed  phase. 
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ni.2  Momentum  equations 

The  equations  of  motion  for  the  fuel  droplet  trajectories  in  the  Cartesian  coordinates  are 


given  below. 

dU. 

pCdA 

2mL 

d  _ 

dt 

AU  (U  -  Ud) 

(m.i) 

dVc 

dt 

pCdA 

2mL 

AU  (V  -  Vd ) 

(HI.2) 

where  the  cross  sectional  area  of  the  fuel  droplet  is  given  by 
A  =  7td2/4 

and  mL  is  the  mass  of  the  droplet.  Ud  and  Vd  are  the  fuel  droplet  velocitiy  components.  The  drag  of 
the  fuel  droplet  is  written  in  terms  of  a  drag  coefficient,  Cd- 

The  relative  velocity  AU  is  given  as 

AU  =  [  (U  -  Ud)2  +  (V  -  Vd)2  ]  ^  (in.3) 

It  is  noted  that  V  (the  y  component  of  gas  phase  mean  velocity)  is  identically  equal  to  zero 
(due  to  the  streamline  orthogonal  coordinate  formulation)  according  to  the  ADD  code  predictions  of 
the  gas  flow  field.  In  the  stochastic  separated  flow  formulation  (see  next  section),  the  contribution 
due  to  the  gas  phase  turbulence  is  considered  in  the  above  equations  of  motion  and  U  and  V  are 
repalced  by  u  and  v  where  (assuming  isotropic  turbulence) 

u  =  U  +  u’ 
and  v  =0  +  u' 

The  fluctuating  component  u'  is  determined  by  random  sampling  of  gas  phase  turbulence  properties 
while  droplet  motion  is  calculated  using  a  random  walk  procedure  (see  next  section) 


The  drag  coefficient  is  determined  from  correlations  reported  by  Dickerson  and  Schuman 

[1965]: 
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Cd  =  27  Red-0-84 
Cd  =  0.271  Red0-217 
Cd  =  2.0 


0  <  Red  <  80 
80  <  Red  <  104 
Red  >  104 


where  the  droplet  Reynolds  number  is  based  on  film  properties  and  the  relative  velocity 


Red  =  Pmd  AU 


The  above  equations  of  motion  constitute  an  initial  value  problem  where  the  initial  values 


Ud  and  Vd  ,  d  and  mL  need  to  be  specified. 


HI.3  Heat  and  mass  transfer  equations  for  droplets 

In  the  case  of  a  (single  component)  droplet  vaporizing  in  a  gas  stream  of  uniform 
temperature,  the  droplet  temperature  history  can  be  described  in  terms  of  a  heat-up  period,  during 
which  the  droplet  temperature  changes,and  an  'equilibrium  vaporization'  period,  during  which  time 
the  droplet  temperature  remains  constant.  The  physical  processes  of  major  significance  during  these 
periods  are  the  transport  of  heat  to  the  droplet  and  the  transport  of  vapor  away  from  the  droplet 
surface..  The  net  difference  between  the  energy  flux  to  and  away  from  the  droplet  accounts  for 
changes  in  the  droplet  temperature. 

The  rigorous  mathematical  treatment  of  the  droplet  vaporization  and  heating  problem 
requires  the  solution  of  time  dependent  partial  diffemtial  equations  for  the  diffusion  of  mass  and 
energy  for  conditions  both  within  the  droplet  and  in  the  surrounding  gas  phase.  Much  progress  has 
been  made  lately  in  the  above  area  where  [Patnaik  ,  1986]  vaporization  from  a  single  droplet  was 
considered.  Since  hundreds  of  classes  of  droplets  are  being  considered  in  a  duct  where  it  is  not  clear 
if  sufficient  information  is  available  for  calculating  the  heat  and  mass  transfer  coefficients  for 
droplets  in  forced  convection,  a  simplified  approach  was  used  in  the  present  research  utilizing  the 
PTRAK  code.  An  existing  semi-empirical  model  [El-Wakil  et  al.,  1954  and  Priem  and  Heidmann, 
1 960]  was  extended  for  use  with  distillate  fuels  (or  pure  compounds)  in  the  PTRAK  code. 

The  following  assumptions  were  made  to  make  the  mathematical  treatment  of  the  droplet 
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heat  and  mass  transfer  tractable: 


(a)  At  each  time  t  during  droplet  vaporization  and  heating,  the  fluid  and  thermodynamic 
conditions  of  the  droplet  and  gas  system  adjust  instantaneously  to  steady  state  conditions  (the 
quasi-steady  assumption).  The  quasi-steady  assumption  is  not  an  equilibrium  assumption. 

(b)  Conditions  within  the  droplet  are  uniform.  This  permits  the  governing  partial  differential 
equations  for  determining  droplet  temperature  and  mass  (or  diameter)  to  be  reduced  to 
extremely  simple  forms  for  quick  solution. 

(c)  Both  the  droplet  and  its  surrounding  are  spherically  symmetric. 

The  geometry  of  the  system  under  consideration  is  shown  in  Figure  17  .  The  droplet  is 
surrounded  by  a  boundary  layer  of  uniform  thickness  and  properties. 


eNtT  —  net  heat  transfer  io  liquid  phasi 
<Ss  -  HEAT  TRANSFER  TO  DROPlET  SURI  ACt 
qA  —  HEAT  REQUIREO  TO  VAPORIZE  FUEL 

<3sn  —  Superheat  rlquireo  to  heat  fuel  vapor  to  ambient  temperature 
4,  —  total  heat  transfer  to  System 

Figure  17.  Geometry  of  droplet  vaporization  model 
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At  the  outer  edge  of  this  boundary  layer,  gas  conditions  correspond  to  those  of  the  flow  field 
calculated  by  the  ADD  code  and  will  be  refertred  to  as  the  "air"  conditions.  The  equation  for  droplet 
vaporization  rate  is  [El  wakil  et  al.,  1954]: 


=  ^Pf.*  a 

where 


Pa 

Pf,S 


ln( 


Pa 

Pa-Pf 


) 


(.111.4) 


(m.5) 


As  is  the  droplet  surface  area,  k  is  the  mass  transfer  coefficient ,  pa  is  the  gas  phase  pressure  and 
Pf>s  is  the  vapor  pressure  at  droplet  surface  at  droplet  temperature  Tl. 


The  mass  transfer  coefficient  k  is  obtained  from  a  correlation  of  the  form  developed  by 
Ranz  and  Marshall  [  1952]  for  droplet  vaprization: 


NUm 


kdTmR 

TmHn 


a  +  bReJ/2  Sc1/3 


(III. 6) 


where  rm  is  the  mass  diffusivity  ,  Mm  is  the  mean  molecular  weight  of  the  air-vapor  mixture 

and  R  is  the  universal  gaS  constant.  The  first  term  on  the  right  hand  side  expression  represents  the 
mass  transfer  due  to  free  convection  while  the  second  term  represents  the  mass  transfer  due  to 
forced  convection.  The  Schmidt  number  is  defined  as: 


Sc=  Prr/(Pmrm) 


The  thermal  coefficients  used  in  these  parameters  are  evaluated  at  a  mean  boundary  layer  film 
condition. 


Mean  film  conditions 

At  any  instant  during  the  evaluation  of  the  vaporization  and  heating  rates,  it  is  assumed 
that  the  appropriate  properties  are  constant  throughout  the  film  surrounding  the  droplet.  The  mean 
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temperature,  Tm,  used  in  the  PTRAK  code  is  the  log  mean  value 


T  =  Ta  ~  Tl  (HI. 7) 

m  ln(T,-TL) 

It  is  also  assumed  that  the  fuel  vapor  and  air  thermal  coefficients  are  obtained  at  a  mean 
fuel  vapor  concentration  (mole  fraction): 

Ym  =  l/2(  Yf  o,  +  Yf  s)  (m.8) 

where  Y  is  the  fuel  vapor  concentration  at  the  droplet  surface. 

The  mole  fraction  is  also  equal  to  the  ratio  of  partial  pressure  to  local  static  pressure  so  that: 

Yf,oo  =  0 

Yf,s  =  Pf,s/Pa 

The  mean  molecular  weight  of  the  mixture  at  the  mean  condition  is: 

Mm  =  YmMf  +(l-Ym)Ma  (m.9) 

and  the  mean  mass  fraction  of  the  fuel  vapor  is 

Cm  =  YmMf/Ma  (III.  10) 

In  addition,  the  five  mean  thermal  properties  that  must  be  evaluated  are  density,  (pm), 
molecular  viscosity  (|im),  thermal  conductivity  (Km),  heat  capacity  (Cpm)  and  the  diffusion 
coefficient,  rm.  From  the  ideal  gas  equation  of  state: 

Pm  =  PaMni/(RTm)  (HI.11) 
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The  mean  molecular  viscosity  (jim)  and  the  mean  thermal  conductivity  are  calculated  as 
mole  weighted  average  quantities.  The  mean  heat  capacity  is  calculated  as 

Cpm  =  On  Cpf  +  (1  •  Cm)  Cpa  (III.  12) 

where  the  heat  capacities  Cpf  and  Cpa  are  evaluated  at  Tm- 


The  diffusion  coefficient  Tm  is  estimated  as  suggested  by  El  Wakil  [1954]  using  the 

kinetic  theory  result  [Hirschfelder  et  al.,  1954].  Force  constants  for  Jet-A  fuel  are  as  given  by 
Anderson  et  al  [1982].  Force  constants  for  water  (as  the  liquid  phase)  and  air  were  obtained  from 
standard  references 


Droplet  heat  transfer 


The  droplet  heating  rate  equation  is 


O 

qs 


WLX 


r  *  Plcpl 


where  the  total  amount  of  heat  arriving  at  the  droplet  surface  is 


(III. 13) 


qs  =  h  As  (Ta  -  Tl)  z/  (ez  -  1)  (III.14) 

and  the  parameter  z  is 

^LCPf  (HI.  15) 

Z  h  As 

w'l  is  the  liquid  fuel  vaporization  rate  and  X  is  the  heat  of  vaporization.  The  droplet  surface  area  As 
is  calculated  from: 


As  =  Jtd^ 

The  above  equation  is  the  product  of  the  forced  convection  heat  transfer  rate  for  a  nonvaporizing 
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system  h  As  (Ta  -  TjJ) ,  and  a  correction  term  for  the  energy  carried  away  from  the  droplet  due  to 
mass  transfer  (i.e.,  "blowing"  of  the  boundary  layer) ,  z  /  (ez  -  1).  This  term  approaches  unity  for 
a  zero  mass  transfer  rate  and  zero  for  an  infinite  mass  transfer  rate  (no  energy  reaching  the  droplet 
surface). 


The  heat  transfer  coefficient  h  is  obtained  from  a  correlation  (similar  to  that  used  for  the 
mass  transfer  coefficient)  developed  by  Ranz  and  Marshall  [1954]: 


Nu  =  ^-=  a’  +  b'Re1/2  Pr1/3  (III. 16) 

h  « 

where  the  Prandtl  number  is  defined  as 
P1  =  ^m^pm^m 

For  reasons  similar  to  those  presented  during  the  discussion  of  the  mass  transfer  coefficient,  the 
coefficients  a'  and  b'  are  assigned  the  values: 


a  -  2.0 
b'  =  0.6 


Some  of  the  heat  arriving  at  the  droplet  surface  is  used  to  vaporize  an  arr,  'int  of  fuel, 

O  O  # 

wl  while  the  remainder  is  used  to  heat  the  droplet.  At  equilibrium,  qs  is  equal  to  wlX  and  the 
droplet  temperature  (the  wet  bulb  temperature)  is  constant  (for  a  pure  fluid). 


At  any  instant,  the  droplet  mass  is  given  by 


mL  =  £  d3  Pl 
Since 


dm 

=  -  w  . 
dt  L 


(111.17) 


(III. 18) 


the  equation  for  the  rate  of  change  of  droplet  diameter  is  given  as 
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(m.i9) 


dt 


nd2PJ 


°  HdLfPk  £5 l 

L  +  6  dTL  dt 


Fuel  properties  required  in  the  model 

Both  liquid  phase  and  vapor  phase  properties  are  required  in  the  droplet  evaporation 
calculations.  The  fuel  vapor  molecular  viscosity  (jif),  the  thermal  conductivity  (Kf)  and  the  heat 
capacity  (Cpf)  are  obtained  (as  function  of  the  mean  film  temperature)  from  polynomial 
expressions. 

For  distillate  fuels,  the  vapor  pressure  (pf>s)  is  evaluated  using  Cox's  procedure  [1923]. 

A  distillate  fuel  is  a  liquid  mixture  consisting  of  two  or  more  compounds  that  may  differ  chemically 
and  physically.  For  a  pure  substance,  such  as  water  or  a  single  hydrocarbon  compound,  the  vapor 
pressure  is  a  function  of  only  the  temperature  of  the  liquid.  For  a  distillate,  it  is  not  practical  from  a 
computational  standpoint  to  keep  track  of  each  component  and  calculate  its  contribution  to  the  vapor 
pressure.  In  Cox's  method  it  is  assumed  that  at  any  instant  during  the  vaporizing  process,  the 
behavior  of  the  distillate  fuel  is  the  same  as  that  of  some  pure  substance;  that  is  its  vapor  pressure  is 
a  function  of  the  liquid  temperature  only.  Different  pure  substances  characterize  the  vaporization 
process  at  each  instant.  Further  details  of  application  of  the  above  procedure  to  the  PTRAK  code 
can  be  found  in  Anderson  et  al.  [1982]. 


The  liquid  density  (p^)  and  heat  capacity  (CpL>  are  calculated  from  polynomial 

expressions  (as  a  function  of  liquid  temperature).  The  liquid  molecular  viscosity  (p^)  is  determined 
from: 


N  =  1 

where  Ajq  are  the  polynomial  coefficients 


(III. 20) 
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Finally,  the  heat  of  vaporization  (X)  is  calculated  using  the  Clausius-Clapeyron  equation  which 
relates  X  to  the  slope  of  the  vapor  pressure  curve: 


d  (  ^  PftS )  -XMf 

dt  rtJ 


(III. 21) 
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m.4  Droplet  class  and  distribution  functions 


For  the  present  calculations  employing  the  PTRAK  code,  the  state  of  a  droplet  at  any 
instant  of  time  is  uniquely  determined  by  its  location  (x,  y),  velocity  (U^  ,  Vj)  in  space,  its 

diameter  (d)  and  its  temperature  (Tl).  In  principle  it  is  is  possible  to  solve  this  set  of  equations  for 
each  of  the  individual  droplets.  However,  the  number  of  droplets  to  be  considered  can  be 
substantially  reduced  by  treating  each  droplet  as  though  it  were  a  cloud  of  'n'  fuel  droplets,  the 
mean  properties  of  which  are  equal  to  the  properties  of  the  individual  droplet. 

A  droplet  class  is  defined  in  the  Lagrangian  sense  by  specifying  the  number  of  droplets  nj 

which  are  within  given  bounds  of  location,  velocity,  droplet  diameter  and  temperature  at  time  =  0. 
In  the  absence  of  droplet  collisions,  the  number  of  droplet  classes  in  phase  space  at  any  given  time 
determines  the  state  of  the  entire  cloud  of  droplets.  The  fuel  flow  rate  is  thus  given  as 

^L  =  1  njirij  (IH.22) 

i  =  o 


Two  distribution  functions  ( the  binomial  and  the  Rosin-Rammler)  were  employed  to 
characterize  the  droplet  sizes  at  the  inlet.  For  the  binomial  distribution  function,  the  mean  particle 
diameter  and  the  particle  variance  (standard  deviation)  need  to  be  specified  along  with  the  number  of 
size  classes.  The  binomial  distribution  function  is  given  as 

f°(1)°  UOL-I)!  (2>1L’  ‘  =0’1L  <In'23) 

where  IL  is  the  number  of  size  classes.  Also 


IL 

X  fp( 0  =  1.0  (III.24) 

>  rr  0 

As  IL  — »  oo  ,  the  binomial  distribution  function  approaches  the  normal  (Gaussian)  distribution 
function. 
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The  Rosin  -  Rammler  distribution  function  is  more  commonly  used  to  characterize  the 
droplet  size  distribution  typically  produced  by  gas  turbine  fuel  injectors.  The  Sauter  mean  diameter 
(SMD)  of  the  spray ,  the  width  factor  m  and  the  number  of  size  classes  have  to  be  specified  for  the 
above  distribution  function.  The  mean  droplet  size  can  be  uniquely  determined  if  the  SMD  and  the 
width  factor  for  the  spray  is  specified.The  width  factor  provides  a  measure  of  the  spread  of  the 
drop  sizes.  The  higher  the  value  of  m,  the  more  uniform  is  the  spray.  If  m  is  infinite,  the  drops  in 
the  spray  are  all  the  same  size.  For  most  fuel  sprays,  the  value  of  m  lies  between  2  and  4.  For 
airblast  atomizers,  Simmons  [1977]  quotes  a  value  of  2.59.  In  the  present  calculations,  the  value  of 
m  was  set  equal  to  3.46. 


III.5  Numerical  Methods 

The  equations  of  motion  for  the  droplets  along  with  the  differential  equations  for  droplet 
diameter  and  temperature  are  solved  as  an  initial  value  problem  using  a  predictor  corrector  method 

which  is  second  order  accurate  in  At.  The  length  of  the  duct  is  divided  into  stations  (usually  finer 
than  the  division  used  in  the  ADD  code)  and  the  initial  guess  for  the  variables  at  station  I  +  1  are  the 
corresponding  values  at  station  I.  The  differential  equations  are  then  evaluated  at  the  mid  point.  The 
equations  are  solved  for  the  unknown  time  (required  for  the  droplets  to  arrive  from  one  station  to 
the  next)  in  an  iterative  manner. 

Depending  on  the  number  of  total  classes  (N)  of  droplets  (product  of  the  droplet  classes 
based  on  initial  diameter,  velocity,  location  and  temperature),  N  sets  of  equations  of  motion  along 
with  the  differential  equations  of  droplet  size  class  temperature  and  droplet  size  class  diameter  are 
solved  in  a  forward  marching  procedure. 
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III.6  Results  and  discussion 


The  fuel  evaporation  effects  in  a  square  duct  were  studied  where  the  fuel  is  injected  from 
the  end  of  an  airblast  atomizer  located  at  the  tip  of  the  splitter  plate  (see  Figure  2).  The  PTRAK  code 
was  employed  (in  conjunction  with  the  ADD  code)  for  this  purpose.Droplet  shattering  and 
coalescence  were  not  considered  in  the  present  study. 

Air  velocities  at  the  exit  plane  of  the  atomizer  are  nonuniform  and  hence  a  turbulent  shear 
layer  forms  downstream  .The  flow  field  in  the  square  duct  is  calculated  first  with  the  modified 

ADD  code  using  a  k  -e  turbulence  model.  Calculations  of  the  nonequilibrium  heating  and 
vaporization  of  the  fuel  droplets  are  then  carried  out  using  the  gas  flow  solution  obtained  first.  In 
this  case  ,  direct  comparisons  of  the  predicted  results  were  not  possible  due  to  the  lack  of 
corresponding  measurements  in  the  parallel  experimental  work.  However,  qualitative  comparisons 
of  the  predictions  and  the  laser  diffraction  measurements  of  the  water  droplets  in  the  tunnel  were  in 
agreement. 

In  this  section,  results  obtained  with  the  unmodified  PTRAK  code  (employing  a 
deterministic  separated  flow  model)  are  presented.  Finite  transport  rates  are  considered  so  that  the 
fuel  droplet  behavior  is  affected  by  the  airflow  behavior.  The  fuel  spray  at  the  exit  of  the  atomizer  is 
characterized  either  by  the  binomial  or  the  Rosin-Rammler  distribution  function.  The  spray  is 
described  by  the  behavior  of  a  number  of  droplet  classes,  each  having  a  different  initial  state 
regarding  the  diameter,  temperature,  velocity  relative  to  the  gas  stream  and  spatial  location. 

The  predictions  of  two  sets  of  calculations  are  reported  below.  In  the  first  set.  Jet  -  A 
(multi-component)  fuel  was  considered  as  the  liquid  phase  and  the  binomial  distribution  function 
was  used  to  characterize  the  spray  at  the  exit  of  the  airblast  atomizer.  In  the  second  set,  a  water 
spray  was  considered  which  was  characterized  by  a  Rosin-Rammler  distribution.  It  is  noted  that  for 
the  experimental  measurements  in  the  parallel  experimental  program,  water  was  atomized  and 
sprayed  in  the  duct. 

For  the  first  case,  uniform  air  velocities  of  60  m/s  (Uafs)  and  30  m/s  (Uaas)  were 
considered  on  the  two  sides  of  the  splitter  plate.  In  order  to  study  the  effects  of  air  flow  field 
characteristics  on  Jet-A  fuel  evaporation,  two  sets  of  inlet  static  pressure  and  stagnation  temperature 
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values  (5  atm. ,  400  K  and  10  atm.  and  750  K)  were  considered  for  the  air  flow.  A  duct  length  of 
2  cm.  was  considered  for  the  computation  of  the  two  dimensional  flow  field.  Insulated  walls  were 
considered  for  the  duct.  For  the  Jet-A  fuel,  the  90%  distillation  temperature  at  1  atm.  was 
considered  to  be  524  K  [Anderson  et  al.,  1982].  A  50  x  40  (  x  x  y)  grid  was  considered  for  the 
computation  of  the  gas  flow  field.  For  the  lower  pressure  -temperature  (5  atm.,  400  K)  gas  flow 
field,  the  Reynolds  number  (based  on  the  width  of  the  duct  and  the  average  inlet  velocity  [Uafs  + 

Uaasl/2  )  was  4.7  x  10^.  The  BaZ  flow  field  predictions  are  similar  to  those  presented  in  an  earlier 
chapter.  The  results  show  the  development  of  a  shear  layer  and  at  larger  values  of  the  downstream 
location,  the  flow  becomes  essentially  self-similar. 

The  Jet-A  fuel  thermophysical  properties  are  given  in  Anderson  et  al.  [1982].  Binomial 
size  distribution  of  the  droplet  classes  were  considered  for  the  calculations.  The  mean  droplet  size  of 
the  spray  was  specified  and  the  standard  deviation  was  set  equal  to  zero.  This,  in  effect  reduced  the 
spray  to  a  monodispersed  spray.  The  fuel  droplets  were  divided  into  ten  classes  (based  upon  the 

injection  angle  and  inlet  velocity  componenets)  between  30°  and  50°  and  -30°  and  -50°  .The  angles 
are  measured  from  the  horizontal  axis,  with  the  anti-clockwise  direction  being  positive.  Each  class 
was  assigned  the  same  initial  temperature.  The  droplets  were  assumed  to  have  an  inlet  velocity  of  30 
m/sec.  The  total  fuel  flow  rate  was  kept  at  0.015  kg/sec  to  ensure  the  low  equivalence  ratios 
necessary  for  the  assumptions  in  the  formulation  of  the  problem. 

The  results  of  these  calculations  can  be  found  in  Farouk  et  al.  [1987]  and  are  not  repeated 
here.  The  purpose  of  these  calculations  were  to  essentially  demonstrate  the  predictive  abilities  of  the 
PTRAK  code  in  the  two  dimensional  duct  situation. 


In  the  second  case  water  was  considered  as  the  dispersed  phase.  For  predictions  of  droplet 
trajectory  and  evaporation  characteristics,  the  gas  flow  field  solution  was  first  obtained  with  Uafs  = 
91.4  m/sec  and  Uaas  =  45  m/sec.  The  inlet  gas  temperature  was  maintained  at  400  K  and  the 
pressure  was  set  equal  to  1  atm.  The  initial  conditions  for  the  gas  phase  computations  were  obtained 
from  measurements  .  The  agreement  of  the  experimental  and  computational  results  of  the  mean 
velocity  and  turbulent  intensity  at  downstream  locations  was  reasonable. 

Water  spray  evporation  was  then  considered  in  the  above  flow  field.  The  spray  was 
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characterized  at  the  inlet  by  the  Rosin-Rammler  distribution  function.  The  initial  SMD  of  the  spray 

was  considered  to  be  91  pm  and  the  spray  was  divided  into  eleven  size  classes.  The  width  factor  of 
the  Rosin-Rammler  distribution  was  set  equal  to  3.46.  To  realistically  simulate  the  conditions  at  the 

exit  of  the  prefilming  atomizer,  six  trajectories  were  considered  (±  1°  to  ±  9°  with  increments  of 

4°)  for  the  spray.  The  eleven  size  classes  were  considered  for  each  trajectory.  Hence  sixty  six  sets 
of  droplet  trajectory,  diameter  and  temperature  equations  were  considered  for  the  solution.  Each 
class  was  assigned  a  small  initial  velocity  (0.1  m/sec)  and  the  initial  droplet  temperatures  were  set 

equal  to  350°  K.The  total  liquid  flow  rate  was  maintained  at  0.015  kg/sec.  The  above  parameters 
for  the  spray  were  similar  to  the  conditions  realized  in  the  parallel  experiments. 

The  calculations  were  carried  out  for  a  duct  length  of  25  cm  (measured  from  the  exit  plane 
of  the  atomizer/splitter  plate).  For  the  length  of  the  computational  domain  and  the  temperatures  of 
the  gas  flow  and  the  incoming  droplets,  the  evaporation  predicted  was  insignificant.  Similar 
observations  were  also  made  in  the  experimental  studies.  Significant  evaporation  of  the  water 
droplets  were  observed  if  the  gas  flow  was  only  at  higher  pressures  and  temperatures.  However, 
those  values  of  pressure  and  temperature  could  not  be  attained  in  the  experimental  facilities 
developed.  The  droplet  velocities  and  temperatures,  were  also  not  measured  in  the  parallel 
experiments. 


Figure  18  shows  the  streamwise  velocities  of  the  largest  (240  |im),  smallest  (20  |im  )  and 

the  intermediate  (105  |im  )  size  class  droplets  injected  at  1°  (towards  the  high  velocity  side)  as  a 
function  of  time.  The  smallest  size  class  attains  the  gas  flow  velocity  within  a  short  time  and  exits 
the  computational  domain  earlier  than  the  larger  size  classes.  The  intermediate  and  the  largest  size 
class  droplets  exhibit  significant  slip  between  their  velocities  and  the  gas  phase  velocity  (Uafs  = 

91.4  for  this  case).  The  temperatures  histories  of  the  above  three  classes  (injected  at  1°  )  are 
shown  in  Figure  19  .  For  all  three  classes,  droplet  temperature  was  found  to  drop  due  to 
evaporation  .  For  the  larger  two  size  classes,  the  temperature  drops  monotonically,  whereas  for  the 
smallest  size  class  droplet  it  is  found  that  due  to  reduced  evaporation,  at  one  point,  the  convective 
heat  transfer  from  the  hot  the  gases  overtakes  the  cooling  due  to  evaporation  .  Similar  observations 
were  made  for  other  size  classes  in  the  other  trajectories  considered.  Due  to  the  lower  velocities  of 

the  gas  phase  flow  encountered  by  trajectories  at  -  1°  through  -  9°,  longer  residence  times  of  the 
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droplets  size  classes  were  encountered. 


The  calculations  were  also  repeated  for  a  different  initial  droplet  temperature  (314°  K) 
while  all  other  parameters  remained  as  above  for  the  gas  and  liquid  phases  [Farouk  et  al. ,  1987b] 


TIME  (MS) 

Figure  18.  Streamwise  droplet  size  class  velocities  for  trajectories  at  1°  angle. 
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19.  Droplet  size  class  temperatures  for  trajectories  at  1°  angle. 


1987b).  The  results  indicated  that  the  amount  of  evaporation  for  the  droplets  was  further  reduced 
and  from  the  very  beginning  all  size  classes  showed  an  increase  of  temperature  . 

Due  to  the  absence  of  the  normal  component  of  the  gas  phase  velocity  (  V  =  0.0)  in  the 
ADD  code  predictions  and  the  deterministic  separated  flow  model  employed  in  the  PTRAK  code, 
the  prediction  of  the  droplet  size  class  trajectories  become  entirely  dependent  on  the  initial  droplet 
velocity  and  initial  trajectory  angle.  Both  of  these  parameters  need  to  be  determined  experimentally 
for  the  prefilming  airblast  atomizer. 


IV.  ANALYSIS  OF  DROPLET  BEHAVIOR -STOCHASTIC  MODEL 
IV.  1  Background  and  general  approach 

Recent  models  for  dilute,  particle-laden  jets  differ  in  details,  but  generally  may  be  divided 
into  two  main  categories  [Shuen,  1984]:  (1)  locally  homogeneous  flow  (LHF)  models,  where 
infinitely  fast  interphase  transport  rates  are  assumed;  (2)  separated  flow  (SF)  models,  where  finite 
interphase  transport  rates  are  considered. 

In  the  Locally  Homogeneous  Flow  (LHF)  models  the  slip  between  the  phases  is 
neglected.  The  main  assumption  of  LHF  models  is  that  interphase  transport  rates  are  fast  in 
comparison  to  the  rate  of  development  of  the  flow  as  a  whole.  This  implies  that  all  phases  have  the 
same  velocity  and  temperature  and  that  phase  equilibrium  is  maintianed  at  each  point  in  the  flow, 
i.e.,the  flow  is  locally  in  thermodynamic  equilibrium.  LHF  models  are  only  correct  when  the 
dispersed  phase  has  infinitely  small  particle  sizes. 

Faeth  and  his  coworkers  studied  evaporating  and  combusting  sprays  [Faeth,  1979,  1980], 
and  condensable  vapor  or  gas  jets  in  liquid  using  LHF  models.  It  was  found  that  while  LHF 
models  provide  a  reasonable  first  estimation  of  the  structure  of  these  processes,  they  generally 
overestimate  the  rate  of  development  of  the  flow.  For  their  experimental  conditions,  a  spray  having 

a  Sauter  mean  diameter  (SMD)  smaller  than  10|i.m  would  be  required  for  quantitative  accuracy  of  the 
LHF  model.  However,  it  was  indicated  that  LHF  models  still  have  some  advantages  [Faeth,  1979], 
namely, 

(1)  they  require  minimal  information  concerning  initial  conditions  of  particle  properties  and 
particle  sizes  -  properties  which  are  often  difficult  to  obtain  for  practical  flows; 

(2)  the  theoretical  model  of  the  flow  is  equivalent  to  that  of  a  single-phase  flow  and  effects  of 
multiple  phases  only  appear  in  the  representation  of  state  relationships  (thermodynamic  properties 
and  molecular  properties),  bypassing  the  difficulties  involved  with  modeling  interactions  between 
the  phases;  and 

(3)  LHF  models  involve  less  empiricism  than  the  SF  models. 


Recent  models  of  dilute  particulate  flows  have  mainly  employed  the  discrete  particle 


formulations.  These  formulations  generally  follow  an  approach  called  "particle- source-in -cell" 
(PSIC)  model  [Crowe,  1977].  This  involves  dividing  the  discrete  phase  into  representative  samples 
whose  motion  and  transport  are  tracked  through  the  flow  field  using  a  Lagrangian  formualtion.  The 
samples  of  the  discrete  phase  are  chosen  to  provide  a  statistically  significant  representation  of  the 
process.  An  Eulerian  formulation  is  generally  employed  to  solve  the  governing  equations  for  the 
continuous  phase,  similar  to  the  solution  of  the  flow  equations  for  LHF 

models.  The  effects  of  interphase  transport  is  considered  by  introducing  appropriate  source  terms 
in  the  governing  equations  for  the  continuous  phase  -  determined  from  the  particle  tracking 
computations.  Most  discrete  particle  models  apply  the  dilute  flow  approximation,  i.e.,  particle 
collisions,  effects  of  nearby  particles  on  particle  transport  rates,  and  the  volume  occupied  by 
particles  are  usually  ignored. 

Deterministic  Separated  Flow  (DSF)  models  adopt  the  discrete  particle  formulation 
neglecting  effects  of  turbulence  on  interphase  transport  properties  and  the  dispersion  of  particles  by 
turbulent  fluctuations.  The  solution  method  employed  in  the  PTRAK  code  falls  under  this  category. 
Neglecting  turbulent  particle  dispersion  is  appropriate  when  characteristic  particle  response  times  are 
large  in  comparison  to  characteristic  times  of  turbulent  fluctuations.  However,  few  practical  flows 
satisfy  this  requirement  Even  when  modified  to  account  for  particle  dispersion,  DSF  models  incur 
errors  when  the  mean  properties  of  the  continuous  phase  are  used  to  evaluate  interphase  transport 
rates  [Faeth,  1983].  Earlier  work  has  shown  that  these  errors  are  particularly  significant  during  the 
evaluation  of  particle  drag. 

In  the  Stochastic  Separated  Flow  (SSF)  models,  effects  of  interphase  slip  and  turbulent 
dispersion  are  considered  using  random  sampling  for  gas  phase  turbulence  properties  in  conjunction 
with  random-walk  computations  for  particle  motion  using  Monte  Carlo  methods  [Gosman  and 
Ionnides,  1984;  Shuen,  1984],  Most  practical  particle-laden  flows  exhibit  properties  which  require 
consideration  of  turbulent  dispersion.  Several  recent  studies  of  turbulent  drop  dispersion  use  SSF 
model  which  have  given  promising  results.  Turbulent  dispersion  of  drops  has  generally  been 
neglected  in  past  separated  flow  models  of  sprays.  Neglecting  dispersion  implies  that  drops  follow 
deterministic  trajectories  prescribed  by  their  initial  condition  at  the  injector  exit  and  mean  gas 
properties  throughout  the  flow  field.  In  locally  homogeneous  flow  models  ,  on  the  other  hand, 
drops  disperse  due  to  turbulence  at  the  same  rate  as  the  gas  phase.  Actual  behavior  lies  between 
these  two  limits,  depending  on  the  inertial  properties  of  the  drops  and  the  turbulence  properties  of 
the  gas  phase.  Techniques  for  treating  turbulent  dispersion  are  rapidly  developing  and  some 
representative  methods  are  discussed  in  the  following. 
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Two  main  methods  have  been  employed  to  compute  particle  dispersion  in  sprays.  The 
earliest  investigations  modeled  turbulent  dispersion  using  a  gradient  diffusion  expression  which 
employs  a  turbulent  gas/particle  exchange  coefficient  Some  method  of  estimating  the  gas/particle 
turbulent  exchange  coefficient  is  required.  Hinze  [1975]  reviews  theoretical  and  experimental 
studies  of  gas/particle  exchange  coefficients.  Existing  continuum  analysis  for  gas/particle  exchange 
coefficients  has  a  limited  range  of  application,  since  the  process  must  be  substantially  simplified  in 
order  to  obtain  closed  form  solutions  [Hinze,  1975].  The  main  difficulty  is  that  dispersion  is 
influenced  by  both  particle  and  turbulence  properties,  adding  a  whole  new  dimension  to  the  problem 
of  estimating  a  turbulent  exchange  coefficient.  Lack  of  theoretical  guidance  on  values  of  gas/particle 
exchange  coefficient  is  compounded  by  difficulties  with  existing  measurements.  These  requirements 
are  rarely  satisfied  in  practice.  Even  in  cases  where  the  experiments  were  properly  analyzed,  it  is 
unlikely  that  the  results  could  be  direcdy  applied  to  a  particular  fuel  spray  since  both  particle  and 
turbulence  properties  must  be  the  same,  because  both  sets  of  characteristics  influence  dispersion 
properties.  Gradient  diffusion  models  of  particle  dispersion  are  not  particularly  convenient  for 
incorporation  in  the  Langrangian  particle  motion  calculatins.  Jurewicz  and  Stock  [1976]  proposes 
circumventing  this  problem  by  representing  the  effects  of  particle  diffusion  as  an  effective  drift 
velocity  which  is  found  from  the  gradient  diffusion  expression  [Jurewicz,  1976].  Unfortunately, 
this  introduces  errors  due  to  numerical  diffusion  which  compromises  one  of  the  advantages  of  the 
Lagrangian  formulation. 

Another  melthod  is  proposed  for  employing  estimated  turbulent  gas/particle  exchange 
coefficients  to  compute  particle  dispersion  [Hotchkiss,  1972,  Dukowicz,  1980].  In  this  case,  the 
turbulent  exchange  coefficient  is  related  to  the  variance  of  the  probability  density  function  of  particle 
position  (assumed  to  be  Gaussian)  at  the  end  of  each  computational  time  step.  Given  the  variance, 
the  distribution  is  randomly  sampled  to  obtain  the  diffusional  increment  of  particle  position.  This 
procedure  does  not  provide  any  real  advantages  over  the  stochastic  methods  [Gosman  and  Ionnides, 
1983],  unless  the  exchange  coefficient  can  be  prescribed  very  simply.  The  method  has  not  been 
applied  to  separated  flow  computations  for  sprays. 

Recent  analysis  of  particle  dispersion  in  sprays  have  employed  stochastic  methods  which 
model  particle  dispersion  directly.  A  turbulence  model  provides  the  instantaneous  properties  of  the 
environment  of  a  particle  during  a  Lagrangian  computation  of  particle  motion.  This  involves 
random  sampling  to  determine  gas  flow  properties  in  addition  to  considering  a  representative 
number  of  particle  trajectories.  Past  computations  with  stochastic  particle  dispersion  models  have 
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used  a  variety  of  turbulence  models.  The  main  features  of  the  stochastic  dispersion  model  of 

Gosman  and  Ioannides  [1983]  are  discussed  in  the  following.  This  method  employs  a  k  -  e 
turbulence  model  and  is  easily  adapted  for  use  in  most  "discrete  droplet  models"  of  sprays. 
Furthermore,  this  model  has  been  employed  to  estimate  the  effect  of  dispersion  in  evaporating  and 
combusting  sprays,  yeilding  resuds  of  some  importance  to  spray  modeling. 

Gosman  and  Ioannides  compute  the  motion  of  a  particle  using  the  Lagrangian  formulation. 
For  particle  dispersion  calculations,  the  instantaneous  gas  velocities  replaces  the  mean  velocities  in 
the  governing  equations.  As  particles  pass  through  the  flow,  they  are  assumed  to  interact  with 
individual  eddies.  Flow  properties  within  each  eddy  are  assumed  to  be  constant.  Each  interaction 
deflects  the  particle  as  dictated  by  the  instantaneous  eddy  velocity.  Therefore,  the  particle  trajectory 
is  determined  similar  to  a  random  walk  computation  until  the  particle  passes  out  of  the  region  under 
consideration.  Overall  particle  behavior  is  obtained  by  averaging  over  a  statistically  significant 
number  of  independent  calculations  . 

IV.2  Mathematical  model 

The  above  stocahstic  model  [  Gosman  and  Ioannides  ,  1983]  for  droplet  dispersion  has 
been  considered  in  the  present  research  for  the  prediction  of  evaparating  sprays  in  a  confined  shear 
layer.  Significant  modifications  to  the  PTRAK  code  has  been  done  to  achieve  this  goal.  The 
instantaneous  gas  velocity  within  an  eddy  (at  a  given  location  )  is  obtained  from  the  mean  gas 
velocity  and  the  turbulent  kinetic  energy  k,  which  are  known  from  the  solution  of  the  governing 
equations  of  the  gas  phase  (using  the  ADD  code).  The  turbulence  is  assumed  to  be  isotropic  with 
fluctuating  components  having  a  Gaussian  distribution  for  this  computation.  The  standard  deviation 
of  the  distribution  is  taken  to  be  (2k/3)^/^.  The  distribution  is  randomly  sampled  when  a  particle 
enters  an  eddy  to  obtain  the  instantaneous  velocities  as 

U  =  u  +  u'  and  V=  0  +u' 

A  particle  is  assumed  to  interact  with  an  eddy  for  a  time  taken  as  the  smaller  of  either  the 
eddy  lifetime,  tg,  or  the  transit  time,  tt,  required  for  the  particle  to  traverse  the  eddy.  The  eddy 
lifetime  and  transit  times  are  determined  by  assuming  that  the  characteristic  size  of  the  randomly 
sampled  eddy  is  the  dissipation  length  scale,  Lg,  given  by 
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(IV.  1) 


Le  =  C^4k3/2/e 


(IV.2) 


av.3) 

(IV  .4) 


IV.  3  Modifications  to  the  PTRAK  code 

The  stochastic  separated  flow  (SSF)  model  involves  finding  trajectories  and  evaporating 
characteristics  of  individual  particles  (  using  a  statistically  significant  number  of  independent 
calculations)  as  they  move  away  from  the  injector  and  encounter  a  random  distribution  of  turbulent 
eddies. 

Key  elements  in  the  SSF  model  are  the  method  used  to  specify  eddy  properties  and  the 
time  of  interaction  of  a  particular  eddy.  An  approach  is  proposed  to  use  direct  simualtion  of  these 
properties  [Shuen,  1984].  Properties  are  assumed  to  be  uniform  within  each  eddy  and  to  change 
randomly  from  one  eddy  to  the  next.  Particle  trajectory  computation  method  is  the  same  as  that 
given  in  the  PTRAK  code,  however,  mean-gas  properties  in  the  particle  equations  are  replaced  by 
the  instantaneous  properties  of  each  eddy. 

The  properties  of  each  eddy  are  found  at  the  start  of  particle-eddy  interaction  by  making  a 
random  selection  from  the  probability  density  function  (PDF)  of  velocity.  Velocity  fluctuaitons  are 
assumed  to  be  isotropic  with  a  Gaussian  PDF  having  a  standard  deviation  of  (2k/3)^  and  mean 
components  of  U  and  zero  in  the  x  and  y  directions  respectively.  Considering  zero  mean 
component  of  velocity  in  the  y  direction  is  the  result  of  using  simplifying  assumptions  in  the  ADD 


The  eddy  lifetime  is  then  estimated  as 

tg  —  Le  /  lu'l 

The  droplet  transit  time  is  determined  from  the  following  equation 
tt  =  -x  In  [1  -  Lg  /  ( I  u  -  Ud  I)  ] 
where  x  is  the  particle  relaxation  time 

x  =  4/3  pfd/(Cd  lu  -  Udl ) 
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code.  The  cummulative  distribution  function  for  each  velocity  component  at  a  given  location  is 
constructed  and  sampled.  This  involves  randomly  selecting  a  number  in  the  range  0  -  1  and 
computing  the  instantaneous  velocity  components  from  the  cummulative  distribution  function. 
This  procedure  provides  a  random  selection  of  eddy  properties  which  satisfies  the  PDF  for  velocity 
fluctuations. 

For  implementing  the  stochastic  separated  flow  model,  we  need  to  randomly  select  the 
instantaneous  velocity,  u,  at  a  point  in  space  ( recall  that  u  =  U  +  u’). 


Since  we  have  assumed  that  the  velocity  distribution  function  is  Gaussian,  the  probabilty  density 
function  of  the  velocity  is  given  as 


p(u)  = — ~=r  exp 
aV27t 


-(u  -  a)2 
2c2 


(IV. 4) 


where  a  is  the  mean  and  a  is  the  standard  distribution. 


If  we  construct  the  cumulative  distribution  function  of  u,  all  the  possible  values  of  u  will 
be  in  the  range  of  0  <  p(u)  <  1. .  A  random  number  generator  will  generate  a  value  between  0  and  1, 
say  A  (0  <  A  <  1),  which  corresponds  to 


u(A) 

J*p(u)  du  =  A 

o 


or 


u(A) 


J: 


o/Irt 


exp 


( u  -  a/ 


2c2 


du 


(IV.5) 


(IV. 6) 


With  the  known  mean  quantity  'a'  and  the  standard  deviation  a,  p(u)  is  uniquely  defined  and  A  is 
given  by  the  random  number  generator.  The  only  unknown  in  the  above  equation  is  the  upper  limit 
u(A)  which  can  be  solved  by  numerical  integration  melthods.  Consequently  we  can  randomly 
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simulate  the  instantaneous  velocity  u  at  any  point  in  space. 


The  equations  of  motion  for  the  fuel  droplet  trajectories  presently  coded  into  the  PTRAK 
computer  program  are  given  in  an  earlier  section.  Uj  and  are  the  mean  droplet  velocities  in  the 
streamwise  and  normal  coordinates  respectively.  The  instantaneous  droplet  velocities  (as  calculated 
in  each  indepenedent  calculation  with  the  stochastic  model)  is  given  by 
ud  =Ud  +  ud' 
and 

vd  =  vd  +  vd' 

where  the  mean  droplet  velocities  are  obtained  by  averaging  a  statistically  significant  number  of 
realizations. 

Summarizing,  the  following  procedure  is  adopted  to  implement  the  stochstic  separated 
flow  model  in  the  ADD  and  the  PTRAK  codes: 


a.  The  solutions  for  the  turbulent  kinetic  energy  (k)  and  the  dissipation  rate  (e)  as  obtained  by 
the  ADD  code  are  stored  in  a  data  file  which  is  subsequently  read  by  the  PTRAK  code. 

b.  Depending  on  the  droplet  size  classlocation  in  the  flow  field,  the  corresponding  values  of 
the  turbulent  kinetic  energy  and  the  length  scale  Le  are  calculated  by  interpolation. 

c.  With  the  mean  gas  velocities  and  the  standard  deviation  known  at  the  droplet  size  class 
location,  the  instantaneous  gas  velocities  are  obtained  by  selecting  a  random  number  via  a  random 
number  generator. 

d.  For  droplet  trajectory  and  evapration  characteristics  calculations,  the  mean  gas  velocities  are 
replaced  by  the  instantaneous  gas  velocities. 

e.  When  the  droplets  arrive  at  the  next  sub-station,  the  time  required  for  the  flight  is  compared 
with  te  and  tt.  If  the  droplet  travel  time  is  larger  than  any  of  tg  and  tt,  step  'c'  is  repeated  for  the  new 
droplet  location.  This  procedure  is  continued  till  the  droplets  exit  the  computational  domain. 

f.  Statistically  significant  number  of  independent  calculations  are  performed.  From  them  , 
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averaged  values  of  U^,  V^,  droplet  diameter  class,  droplet  temperature  etc.  are  computed  along  the 
length  of  the  duct. 


67 


IV.4  Results  and  discussion 


The  modified  PTRAK  code  along  with  the  gas  flow  predictions  by  the  ADD  code  can  be 
employed  to  predict  the  spray  evaporation  characteristics  in  the  turbulent  duct  flow.  The  effect  of 
the  gas  phase  turbulence  on  droplet  dispersion  is  considered  via  a  stochastic  separated  flow  model. 
The  amount  of  droplet  evaporation  in  the  parallel  experiments  was  found  to  be  insignificant  for  the 
gas  flow  pressure  and  temperatures  in  the  duct  and  the  water  temperature  at  the  atomizer.  Since  the 
predictions  for  the  stochastic  separated  flow  model  were  obtained  for  conditions  similar  to  those  in 
the  experiments,  the  droplet  evaporation  within  the  25  cm  of  the  duct  length  considered  was  small. 
The  dispersion  effects  of  the  droplet  size  classes  were  only  of  significance.  The  changes  due  to  gas 
phase  turbulence  on  the  droplet  temperature  history  and  the  evaporation  rate  was  minimal,  for  the 
operating  conditions  studied. 

A  wake  flow  field  was  considered  for  gas  phase  computations  where  Uafs  and  Uaas  were 
both  set  equal  to  91.4  m/sec.  The  inlet  values  for  the  mean  axial  velocity  was  obtained  from 
experimental  measurements  at  x/D  =  0.03.  The  pressure  and  temperature  of  the  gas  flow  field  at  the 
inlet  was  uniform  and  was  set  to  1  atm.  and  400  K  respectively.  Due  to  the  presence  of  the  splitter 
plate,  a  turbulent  wake  flow  is  created  downstream.  The  mean  axial  flow  solution,  U(x, y)  along 

with  the  turbulent  kinetic  energy  and  dissipation  rate  solutions,  k(x,y)  and  e(x,y)  were  stored  and 
subsequently  read  by  the  PTRAK  code.  Water  was  considered  as  the  dispersed  phase  and  a 
Rosin-Rammler  distribution  function  was  employed  to  characterize  the  spray  at  the  tip  of  the  airblast 

atomizer.  The  spray  was  divided  into  eleven  size  classes  with  SMD  =  120pm  and  the  width  factor 
of  3.46.  In  order  to  investigate  the  dispersion  characteristics  of  the  spray,  only  one  trajectory  of  the 

spray  was  considered  with  an  injection  angle  of  1  °.  Each  class  was  assigned  the  same  initial 
temperature  of  314  K.  Though  only  one  trajectory  was  considered,  the  fuel  flow  rate  was  still 
maintained  at  to  0.015kg/sec.  For  the  initial  velocity  of  the  droplet  size  classes,  several  cases  were 
investigated.  It  was  observed  that  the  dispersion  effects  due  to  gas  phase  turbulence  diminshes  with 

increasing  inlet  velocities  of  the  droplets.  Figure  20  shows  the  trajectories  of  the  smallest  (24  pm) 

the  largest  (340  pm)  and  the  intermediate  size  (210  pm)  classes  using  the  unmodified  PTRAK  code. 
The  initial  droplet  velocity  was  set  equal  to  10  m/sec.  No  effect  of  the  gas  phase  turbulence  is 
accounted  for  in  this  result.  The  smallest  size  class  exits  the  computational  domain  first  and  is 
strongly  influenced  by  the  parallel  (axial )  gas  flow  field.  The  larger  size  classes  (due  to  the  initial 
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injection  angle  and  larger  inertia)  continue  to  move  obliquely  in  the  parallel  gas  flow  field. 

The  modified  PTRAK  code  (with  a  stochastic  separated  flow  model )  was  then  applied  to 
calculate  the  trajectories  and  evaporation  characteristics  of  the  above  size  classes  with  identical 
initial  conditions.  A  statistically  significant  number  of  independent  computer  runs  were  made  where 
the  instantaneous  gas  flow  field  at  any  location  (u  =  U  +  u’  and  v  =  0  +  u')  was  simulated  by 
procedures  described  earlier.  The  stored  values  of  the  turbulent  kinetic  energy  was  used  to  compute 
the  vallues  of  u’  at  any  droplet  location.  The  stored  values  of  the  the  dissipation  rate  of  turbulence 
was  used  to  compute  the  eddy  lifetime  tg  and  the  droplet  transit  time  t{  as  described  earlier.  Figures 
21,  22  and  23  show  the  droplet  size  class  trajectories  for  three  independent  calculaions  performed 
by  the  modified  PTRAK  code.  Comparing  these  results  with  Fig.  20  (obtained  by  using  the 
unmodified  PTRAK  code),  the  efects  of  the  gas  phase  turbulence  on  the  droplet  trajectories  is 
evident.  As  expected,  the  effects  are  stronger  for  the  smallest  size  droplets.  Due  to  the  stochastic 
approach  used,  the  results  in  Figs.  21,  22  and  23  vary  in  details,  but  qualitatively  exhibit  the  same 
behavior.  A  'statistically  significant*  number  of  these  computations  (realizations)  were  obtained  and 
then  the  solutions  were  averaged  .  Figure  24  shows  the  predicted  trajectories  of  the  droplet  classes 
after  the  averaging  process. 
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Figure  20.  Droplet  size  class  trajectories  obtained  with  the  unmodified  PTRAK  code. 
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Figure  21.  A  realization  of  the  droplet  size  class  trajectories  obtained  with  the  modified  PTRAK 
code  . 
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Figure  22.  A  realization  of  the  droplet  size  class  trajectories  obtained  with  the  modified  PTRAK 
code . 
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Figure  23.  A  realization  of  the  droplet  size  class  trajectories  obtained  with  the  modified  PTRAK 
code  . 
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Figure  24.  Droplet  size  class  trajectories  (averaged)  obtained  with  the  modified  PTRAK  code  . 
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Comfwing  Figures  20  and  24,  it  seen  that  at  the  exit  of  the  computational  domair,  the 
smalles  size  class  droplet  is  deflected  (in  the  nor.mal  direction)  the  highest  amont  due  to  the  gas 
phase  turbulence  (about  135%  more  than  he  predictions  obtained  with  the  deterministic  separated 
flow  model).  For  the  intermediate  and  the  largest  size  classes  the  effects  were  less  pronounced,  but 
significant  (80%  and  45%  more  deflection  in  the  normal  direction).  Similar  behavior  of  the  spray 
was  observed  for  the  spray  in  the  experimental  studies.  Droplet  size  class  temperature  and 
evaporation  predictions  with  the  stochastic  separated  flow  model  did  not  show  significant 
variations  with  those  obtained  with  the  unmodified  PTRAK  code  for  the  flow  and  initial  droplet  size 
conditions  considered.  These  are  primarily  due  to  the  low  temperature  difference  between  the  gas 
phase  flow  and  the  initial  liquid  temperatures.  In  a  gas  turbine  air-fuel  passage,  the  air  is  at  higher 
pressure  which  enhances  fuel  evaporation.  The  effects  of  droplet  dispersion  can  be  studied  more 
effectively  with  the  modified  PTRAK  code  if  the  appropriate  gas  flow  field  conditions  are  used  in 
the  simulations. 
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V. 


SUMMARY  AND  CONCLUSIONS 


A  computational  study  was  undertaken  to  study  the  evaporation  characteristics  of  a  liquid 
spray  in  a  confined  turbulent  shear  flow.  The  computations  modeled  parallel  experiments  carried 
out  under  ARO  sponsorship  (Grant  No.  DAAG29-84-K-0165).  The  experiments  were  carried  out 
in  a  vertically  down-flowing  wind  tunnel  where  a  splitter  plate  separates  two  air  streams  whose 
average  velocities  can  be  controlled  independently.  Depending  upon  the  values  of  the  avearge 
velocities  on  each  side  of  the  splitter  plate,  a  turbulent  shear  layer  or  wake  flow  develops  at  the  end 
of  the  splitter  plate  in  the  duct..  A  flat  prefilming  airblast  atomizer  is  located  at  the  tip  of  the  splitter 
plate.  The  trajectories  and  the  evaporation  characteristics  of  the  water  spray  formed  at  the  atomizer 
are  affected  by  the  turbulent  gas  phase  flow  field.  The  droplet  formation  process  is  also  affected  by 
the  flow  characteristics  on  the  two  sides  of  the  splitter  plate. 

The  computer  codes  ADD  and  PTRAK  [Anderson  et  al.,  1982]  developed  for  the 
prediction  of  flow  field  and  fuel  evaporation  characteristics  in  prevaporizing-premixing  fuel  air 
passages  (of  gas  turbine  combustors)  were  chosen  for  carrying  out  the  computations.  This  was 
perhaps  a  poor  choice  due  to  the  complexities  of  the  turbulent  flow  field  in  the  experimental 
configuration  and  the  spray  formation  process  in  the  airblast  atomizer.  The  flow  situations 
considered  in  the  present  research  program  are  considerably  more  complicated  than  the 
axisymmetric  parabolic  duct  flow  situations  that  are  suitable  for  the  ADD  and  the  PTRAK  codes.  In 
fact,  the  flow  field  within  a  gas  turbine  combustor  (dominated  by  strong  shear  layers)  are  inherently 
more  complex  than  those  found  in  axisymmetric  ducts.  Extensive  modifications  to  the  ADD  code 
had  to  be  performed  to  obtain  predictions  for  the  two  dimensional  turbulent  shear  flows  (at  the  end 
of  a  splitter  plate  located  along  the  axis  of  the  duct)  which  agreed  reasonably  well  with  the 
measurements. 

The  normal  component  of  the  mean  gas  velocity  is  not  solved  by  the  ADD  code  (due  to  the 
orthogonal  streamline  coordinate  system  employed).  Neither  can  it  handle  significant  flow  reversal. 
These  introduce  errors  in  calculating  flow  fields  at  the  exit  of  the  splitter  plate.  Due  to  the  parabolic 
formulation,  it  was  also  not  possible  to  include  the  splitter  plate  in  the  computational  domain. 
Recent  studies  have  shown  that  considering  the  splitter  palte  within  the  computational  domain 
significantly  enhances  the  accuracy  of  the  shear  layer  predictions  [Grinstein  et  al.,  1986].  Current 
efforts  are  underway  at  Drexel  University  to  compute  the  flow  field  in  the  experimental  tunnel  using 

an  elliptic  formulation  of  the  governing  equations.  Efforts  in  tuning  the  parameters  of  the  k  -e 
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model  for  the  gas  phase  turbulence  produced  some  improvements  in  the  predictions. 

The  ADD  and  PTRAK  codes  suffer  from  another  major  limitation  in  the  formulation,  i.e. 
there  is  only  one  way  coupling  between  the  gas  phase  and  the  dispersed  phase  computations.  The 
effect  of  droplet  evaporation  is  not  considered  in  the  gas  flow  computations.  This  assumption  is 
reasonable  for  computations  in  a  long  prevaporizing-premixing  duct.  However,  for  computations 
close  to  an  atomizer  in  a  combustion  chamber,  this  can  introduce  significant  errors. 

The  computations  in  the  PTRAK  code  are  suited  for  droplets  generated  by  a  pressure 
atomizer  where  the  gas  phase  flow  field  has  little  influence  in  the  spray  formation  process.  For  the 
airblast  atomizer  (as  used  in  the  experimental  set  up)  the  spray  formation  process  is  found  to  be 
strongly  influenced  by  the  gas  flow  field,  in  particular,  by  the  shear  layer  strength.  Hence  for 
realistic  prediction  of  the  trajectories  and  the  evaporation  characteristics,  the  droplet  size  class 
properties  need  to  be  characterized  at  the  inlet  which  are  in  effect  dictated  by  the  gas  flow  field  .  For 
the  PTRAK  code,  the  droplet  size  classes  at  the  inlet  need  to  be  characterized  by  their  velocities  and 
the  injection  angles.  These  parameters  are  more  easily  found  for  a  pressure  atomizer  than  an  airblast 
atomizer  where  the  spray  formation  process  is  independent  of  the  gas  flow  velocities.  Future 
research  efforts  should  be  directed  towards  estimating  the  initial  conditions  for  the  droplet  size 
classes  formed  in  an  airblast  atomizer  for  different  flow  conditions. 

The  modifications  made  to  the  PTRAK  code  has,  however,  improved  its  predictive 
capabilities  significantly.  A  major  deficiency  of  the  ADD  code  is  that  the  normal  component  of  the 
mean  gas  phase  velocity  is  set  to  zero  everywhere  in  the  flow  domain.  With  the  addition  of  the 
stochastic  separated  flow  model,  the  effects  of  gas  phase  turbulence  is  considered  in  the  droplet  size 
class  trajectory  and  evaporation  characteristics  calculations.  This  in  effect  introduces  a  non-zero 
normal  component  of  the  gas  phase  velocity  in  the  droplet  trajectory  calculations.  Computations  are 
being  carried  out  at  Drexel  to  predict  the  spray  evaporation  charactersitics  in  a 
prevaporizing-premixing  duct  with  the  modified  PTRAK  code. 
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